Merge branch 'develop' into collected-small-fixes

cmake/CMakeLists.txt

@@ -440,7 +440,7 @@ if(BUILD_OMP)
   target_link_libraries(lmp PRIVATE OpenMP::OpenMP_CXX)
 endif()
 
-if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE)
+if(PKG_MSCG OR PKG_ATC OR PKG_AWPMD OR PKG_ML-QUIP OR PKG_ML-POD OR PKG_ELECTRODE OR BUILD_TOOLS)
   enable_language(C)
   if (NOT USE_INTERNAL_LINALG)
     find_package(LAPACK)

cmake/Modules/Tools.cmake

@@ -33,6 +33,8 @@ if(BUILD_TOOLS)
   endif()
   install(TARGETS msi2lmp DESTINATION ${CMAKE_INSTALL_BINDIR})
   install(FILES ${LAMMPS_DOC_DIR}/msi2lmp.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1)
+
+  add_subdirectory(${LAMMPS_TOOLS_DIR}/phonon ${CMAKE_BINARY_DIR}/phana_build)
 endif()
 
 if(BUILD_LAMMPS_SHELL)

doc/src/Commands_compute.rst

@@ -46,6 +46,7 @@ KOKKOS, o = OPENMP, t = OPT.
    * :doc:`com/chunk <compute_com_chunk>`
    * :doc:`contact/atom <compute_contact_atom>`
    * :doc:`coord/atom (k) <compute_coord_atom>`
+   * :doc:`count/type <compute_count_type>`
    * :doc:`damage/atom <compute_damage_atom>`
    * :doc:`dihedral <compute_dihedral>`
    * :doc:`dihedral/local <compute_dihedral_local>`

doc/src/Howto_broken_bonds.rst

@@ -1,48 +1,56 @@
 Broken Bonds
 ============
 
-Typically, bond interactions persist for the duration of a simulation in
+Typically, molecular bond interactions persist for the duration of a
-LAMMPS. However, there are some exceptions that allow for bonds to
+simulation in LAMMPS.  However, some commands break bonds dynamically,
-break, including the :doc:`quartic bond style <bond_quartic>` and the
+including the following:
-bond styles in the :doc:`BPM package <Howto_bpm>` which contains the
-:doc:`bpm/spring <bond_bpm_spring>` and :doc:`bpm/rotational
-<bond_bpm_rotational>` bond styles. In these cases, a bond can be broken
-if it is stretched beyond a user-defined threshold.  LAMMPS accomplishes
-this by setting the bond type to 0, such that the bond force is no
-longer computed.
 
-Users are normally able to weight the contribution of pair forces to atoms
+* :doc:`bond_style quartic <bond_quartic>`
-that are bonded using the :doc:`special_bonds command <special_bonds>`.
+* :doc:`fix bond/break <fix_bond_break>`
-When bonds break, this is not always the case. For the quartic bond style,
+* :doc:`fix bond/react <fix_bond_react>`
-pair forces are always turned off between bonded particles. LAMMPS does
+* :doc:`BPM package <Howto_bpm>` bond styles
-this via a computational sleight-of-hand. It subtracts the pairwise
-interaction as part of the bond computation. When the bond breaks, the
-subtraction stops.  For this to work, the pairwise interaction must always
-be computed by the :doc:`pair_style <pair_style>` command, whether the bond
-is broken or not.  This means that :doc:`special_bonds <special_bonds>` must
-be set to 1,1,1. After the bond breaks, the pairwise interaction between the
-two atoms is turned on, since they are no longer bonded.
 
-In the BPM package, one can either turn off all pair interactions between
+A bond can break if it is stretched beyond a user-defined threshold or
-bonded particles or leave them on, overlaying pair forces on top of bond
+more generally if other criteria are met.
-forces. To remove pair forces, the special bond list is dynamically
-updated. More details can be found on the :doc:`Howto BPM <Howto_bpm>`
-page.
 
-Bonds can also be broken by fixes which change bond topology, including
+For the quartic bond style, when a bond is broken its bond type is set
-:doc:`fix bond/break <fix_bond_break>` and
+to 0 to effectively break it and pairwise forces between the two atoms
-:doc:`fix bond/react <fix_bond_react>`. These fixes will automatically
+in the broken bond are "turned on".  Angles, dihedrals, etc cannot be
-trigger a rebuild of the neighbor list and update special bond data structures
+defined for a system when :doc:`bond_style quartic <bond_quartic>` is
-when bonds are broken.
+used.
 
-Note that when bonds are dumped to a file via the :doc:`dump local <dump>` command, bonds with type 0 are not included.  The
+Similarly, bond styles in the BPM package are also incompatible with
-:doc:`delete_bonds <delete_bonds>` command can also be used to query the
+angles, dihedrals, etc. and when a bond breaks its type is set to zero.
-status of broken bonds or permanently delete them, e.g.:
+However, in the BPM package one can either turn off all pair interactions
+between bonded particles or leave them on, overlaying pair forces on
+top of bond forces. To remove pair forces, the special bond list is
+dynamically updated.  More details can be found on the :doc:`Howto BPM
+<Howto_bpm>` page.
+
+The :doc:`fix bond/break <fix_bond_break>` and :doc:`fix bond/react
+<fix_bond_react>` commands allow breaking of bonds within a molecular
+topology with may also define angles, dihedrals, etc.  These commands
+update internal topology data structures to remove broken bonds, as
+well as the appropriate angle, dihederal, etc interactions which
+include the bond.  They also trigger a rebuild of the neighbor list
+when this occurs, to turn on the appropriate pairwise forces.
+
+Note that when bonds are dumped to a file via the :doc:`dump local
+<dump>` command, bonds with type 0 are not included.
+
+The :doc:`delete_bonds <delete_bonds>` command can be used to query
+the status of broken bonds with type = 0 or permanently delete them,
+e.g.:
 
 .. code-block:: LAMMPS
 
    delete_bonds all stats
    delete_bonds all bond 0 remove
 
-The compute :doc:`nbond/atom <compute_nbond_atom>` can also be used
+The compute :doc:`count/type <compute_count_type>` command tallies the
-to tally the current number of bonds per atom, excluding broken bonds.
+current number of bonds (or angles, etc) for each bond (angle, etc)
+type.  It also tallies broken bonds with type = 0.
+
+The compute :doc:`nbond/atom <compute_nbond_atom>` command tallies the
+current number of bonds each atom is part of, excluding broken bonds
+with type = 0.

doc/src/Tools.rst

@@ -883,9 +883,9 @@ dependencies and redirects the download to the local cache.
 phonon tool
 ------------------------
 
-The phonon subdirectory contains a post-processing tool useful for
+The phonon subdirectory contains a post-processing tool, *phana*, useful
-analyzing the output of the :doc:`fix phonon <fix_phonon>` command in
+for analyzing the output of the :doc:`fix phonon <fix_phonon>` command
-the PHONON package.
+in the PHONON package.
 
 See the README file for instruction on building the tool and what
 library it needs.  And see the examples/PACKAGES/phonon directory

doc/src/bond_style.rst

@@ -32,13 +32,13 @@ Set the formula(s) LAMMPS uses to compute bond interactions between
 pairs of atoms.  In LAMMPS, a bond differs from a pairwise
 interaction, which are set via the :doc:`pair_style <pair_style>`
 command.  Bonds are defined between specified pairs of atoms and
-remain in force for the duration of the simulation (unless the bond
+remain in force for the duration of the simulation (unless new bonds
-breaks which is possible in some bond potentials).  The list of bonded
+are created or existing bonds break, which is possible in some fixes
-atoms is read in by a :doc:`read_data <read_data>` or
+and bond potentials).  The list of bonded atoms is read in by a
-:doc:`read_restart <read_restart>` command from a data or restart file.
+:doc:`read_data <read_data>` or :doc:`read_restart <read_restart>`
-By contrast, pair potentials are typically defined between all pairs
+command from a data or restart file.  By contrast, pair potentials are
-of atoms within a cutoff distance and the set of active interactions
+typically defined between all pairs of atoms within a cutoff distance
-changes over time.
+and the set of active interactions changes over time.
 
 Hybrid models where bonds are computed using different bond potentials
 can be setup using the *hybrid* bond style.

doc/src/compute.rst

@@ -200,6 +200,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
 * :doc:`com/chunk <compute_com_chunk>` - center of mass for each chunk
 * :doc:`contact/atom <compute_contact_atom>` - contact count for each spherical particle
 * :doc:`coord/atom <compute_coord_atom>` - coordination number for each atom
+* :doc:`count/type <compute_count_type>` - count of atoms or bonds by type
 * :doc:`damage/atom <compute_damage_atom>` - Peridynamic damage for each atom
 * :doc:`dihedral <compute_dihedral>` - energy of each dihedral sub-style
 * :doc:`dihedral/local <compute_dihedral_local>` - angle of each dihedral

doc/src/compute_count_type.rst

@@ -0,0 +1,130 @@
+.. index:: compute count/type
+
+compute count/type command
+====================
+
+Syntax
+""""""
+
+.. code-block:: LAMMPS
+
+   compute ID group-ID count/type mode
+
+* ID, group-ID are documented in :doc:`compute <compute>` command
+* count/type = style name of this compute command
+* mode = {atom} or {bond} or {angle} or {dihedral} or {improper}
+
+Examples
+""""""""
+
+.. code-block:: LAMMPS
+
+   compute 1 all count/type atom
+   compute 1 flowmols count/type bond
+
+Description
+"""""""""""
+
+.. versionadded:: TBD
+
+Define a computation that counts the current number of atoms for each
+atom type.  Or the number of bonds (angles, dihedrals, impropers) for
+each bond (angle, dihedral, improper) type.
+
+The former can be useful if atoms are added to or deleted from the
+system in random ways, e.g. via the :doc:`fix deposit <fix_deposit>`,
+:doc:`fix pour <fix_pour>`, or :doc:`fix evaporate <fix_evaporate>`
+commands.  The latter can be useful in reactive simulations where
+molecular bonds are broken or created, as well as angles, dihedrals,
+impropers.
+
+Note that for this command, bonds (angles, etc) are the topological
+kind enumerated in a data file, initially read by the :doc:`read_data
+<read_data>` command or defined by the :doc:`molecule <molecule>`
+command.  They do not refer to implicit bonds defined on-the-fly by
+bond-order or reactive pair styles based on the current conformation
+of small clusters of atoms.
+
+These commands can turn off topological bonds (angles, etc) by setting
+their bond (angle, etc) types to negative values.  This command
+includes the turned-off bonds (angles, etc) in the count for each
+type:
+
+* :doc:`fix shake <fix_shake>`
+* :doc:`delete_bonds <delete_bonds>`
+
+These commands can create and/or break topological bonds (angles,
+etc).  In the case of breaking, they remove the bond (angle, etc) from
+the system, so that they no longer exist (:doc:`bond_style quartic
+<bond_quartic>` and :doc:`BPM bond styles <Howto_bpm>` are exceptions,
+see the discussion below).  Thus they are not included in the counts
+for each type:
+
+* :doc:`delete_bonds remove <delete_bonds>`
+* :doc:`bond_style quartic <bond_quartic>`
+* :doc:`fix bond/react <fix_bond_react>`
+* :doc:`fix bond/create <fix_bond_create>`
+* :doc:`fix bond/break <fix_bond_break>`
+* :doc:`BPM package <Howto_bpm>` bond styles
+
+----------
+
+If the {mode} setting is {atom} then the count of atoms for each atom
+type is tallied.  Only atoms in the specified group are counted.
+
+If the {mode} setting is {bond} then the count of bonds for each bond
+type is tallied.  Only bonds with both atoms in the specified group
+are counted.
+
+For {mode} = {bond}, broken bonds with a bond type of zero are also
+counted.  The :doc:`bond_style quartic <bond_quartic>` and :doc:`BPM
+bond styles <Howto_bpm>` break bonds by doing this.  See the :doc:`
+Howto broken bonds <Howto_broken_bonds>` doc page for more details.
+Note that the group setting is ignored for broken bonds; all broken
+bonds in the system are counted.
+
+If the {mode} setting is {angle} then the count of angles for each
+angle type is tallied.  Only angles with all 3 atoms in the specified
+group are counted.
+
+If the {mode} setting is {dihedral} then the count of dihedrals for
+each dihedral type is tallied.  Only dihedrals with all 4 atoms in the
+specified group are counted.
+
+If the {mode} setting is {improper} then the count of impropers for
+each improper type is tallied.  Only impropers with all 4 atoms in the
+specified group are counted.
+
+----------
+
+Output info
+"""""""""""
+
+This compute calculates a global vector of counts.  If the mode is
+{atom} or {bond} or {angle} or {dihedral} or {improper}, then the
+vector length is the number of atom types or bond types or angle types
+or dihedral types or improper types, respectively.
+
+If the mode is {bond} this compute also calculates a global scalar
+which is the number of broken bonds with type = 0, as explained above.
+
+These values can be used by any command that uses global scalar or
+vector values from a compute as input.  See the :doc:`Howto output
+<Howto_output>` page for an overview of LAMMPS output options.
+
+The scalar and vector values calculated by this compute are "extensive".
+
+Restrictions
+""""""""""""
+
+none
+
+Related commands
+""""""""""""""""
+
+none
+
+Default
+"""""""
+
+none

src/KOKKOS/fix_langevin_kokkos.cpp

@@ -44,6 +44,7 @@ FixLangevinKokkos<DeviceType>::FixLangevinKokkos(LAMMPS *lmp, int narg, char **a
   FixLangevin(lmp, narg, arg),rand_pool(seed + comm->me)
 {
   kokkosable = 1;
+  fuse_integrate_flag = 1;
   sort_device = 1;
   atomKK = (AtomKokkos *) atom;
   int ntypes = atomKK->ntypes;
@@ -170,6 +171,14 @@ void FixLangevinKokkos<DeviceType>::initial_integrate_item(int i) const
 
 /* ---------------------------------------------------------------------- */
 
+template<class DeviceType>
+void FixLangevinKokkos<DeviceType>::fused_integrate(int vflag)
+{
+  initial_integrate(vflag);
+}
+
+/* ---------------------------------------------------------------------- */
+
 template<class DeviceType>
 void FixLangevinKokkos<DeviceType>::post_force(int /*vflag*/)
 {

src/KOKKOS/fix_langevin_kokkos.h

@@ -70,6 +70,7 @@ namespace LAMMPS_NS {
     void cleanup_copy();
     void init() override;
     void initial_integrate(int) override;
+    void fused_integrate(int) override;
     void post_force(int) override;
     void reset_dt() override;
     void grow_arrays(int) override;

src/KOKKOS/fix_nve_kokkos.cpp

@@ -29,6 +29,7 @@ FixNVEKokkos<DeviceType>::FixNVEKokkos(LAMMPS *lmp, int narg, char **arg) :
   FixNVE(lmp, narg, arg)
 {
   kokkosable = 1;
+  fuse_integrate_flag = 1;
   atomKK = (AtomKokkos *) atom;
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
 
@@ -159,6 +160,66 @@ void FixNVEKokkos<DeviceType>::final_integrate_rmass_item(int i) const
   }
 }
 
+/* ----------------------------------------------------------------------
+   allow for both per-type and per-atom mass
+------------------------------------------------------------------------- */
+
+template<class DeviceType>
+void FixNVEKokkos<DeviceType>::fused_integrate(int /*vflag*/)
+{
+  atomKK->sync(execution_space,datamask_read);
+
+  x = atomKK->k_x.view<DeviceType>();
+  v = atomKK->k_v.view<DeviceType>();
+  f = atomKK->k_f.view<DeviceType>();
+  rmass = atomKK->k_rmass.view<DeviceType>();
+  mass = atomKK->k_mass.view<DeviceType>();
+  type = atomKK->k_type.view<DeviceType>();
+  mask = atomKK->k_mask.view<DeviceType>();
+  int nlocal = atomKK->nlocal;
+  if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst;
+
+  if (rmass.data()) {
+    FixNVEKokkosFusedIntegrateFunctor<DeviceType,1> functor(this);
+    Kokkos::parallel_for(nlocal,functor);
+  } else {
+    FixNVEKokkosFusedIntegrateFunctor<DeviceType,0> functor(this);
+    Kokkos::parallel_for(nlocal,functor);
+  }
+
+  atomKK->modified(execution_space,datamask_modify);
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void FixNVEKokkos<DeviceType>::fused_integrate_item(int i) const
+{
+  if (mask[i] & groupbit) {
+    const double dtfm = 2.0 * dtf / mass[type[i]];
+    v(i,0) += dtfm * f(i,0);
+    v(i,1) += dtfm * f(i,1);
+    v(i,2) += dtfm * f(i,2);
+    x(i,0) += dtv * v(i,0);
+    x(i,1) += dtv * v(i,1);
+    x(i,2) += dtv * v(i,2);
+  }
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void FixNVEKokkos<DeviceType>::fused_integrate_rmass_item(int i) const
+{
+  if (mask[i] & groupbit) {
+    const double dtfm = 2.0 * dtf / rmass[i];
+    v(i,0) += dtfm * f(i,0);
+    v(i,1) += dtfm * f(i,1);
+    v(i,2) += dtfm * f(i,2);
+    x(i,0) += dtv * v(i,0);
+    x(i,1) += dtv * v(i,1);
+    x(i,2) += dtv * v(i,2);
+  }
+}
+
 /* ---------------------------------------------------------------------- */
 
 template<class DeviceType>
@@ -168,6 +229,8 @@ void FixNVEKokkos<DeviceType>::cleanup_copy()
   vatom = nullptr;
 }
 
+/* ---------------------------------------------------------------------- */
+
 namespace LAMMPS_NS {
 template class FixNVEKokkos<LMPDeviceType>;
 #ifdef LMP_KOKKOS_GPU

src/KOKKOS/fix_nve_kokkos.h

@@ -46,6 +46,7 @@ class FixNVEKokkos : public FixNVE {
   void init() override;
   void initial_integrate(int) override;
   void final_integrate() override;
+  void fused_integrate(int) override;
 
   KOKKOS_INLINE_FUNCTION
   void initial_integrate_item(int) const;
@@ -55,6 +56,10 @@ class FixNVEKokkos : public FixNVE {
   void final_integrate_item(int) const;
   KOKKOS_INLINE_FUNCTION
   void final_integrate_rmass_item(int) const;
+  KOKKOS_INLINE_FUNCTION
+  void fused_integrate_item(int) const;
+  KOKKOS_INLINE_FUNCTION
+  void fused_integrate_rmass_item(int) const;
 
  private:
 
@@ -96,6 +101,22 @@ struct FixNVEKokkosFinalIntegrateFunctor  {
   }
 };
 
+template <class DeviceType, int RMass>
+struct FixNVEKokkosFusedIntegrateFunctor  {
+  typedef DeviceType  device_type ;
+  FixNVEKokkos<DeviceType> c;
+
+  FixNVEKokkosFusedIntegrateFunctor(FixNVEKokkos<DeviceType>* c_ptr):
+  c(*c_ptr) {c.cleanup_copy();};
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int i) const {
+    if (RMass)
+      c.fused_integrate_rmass_item(i);
+    else
+      c.fused_integrate_item(i);
+  }
+};
+
 }
 
 #endif

src/KOKKOS/fix_nve_sphere_kokkos.cpp

@@ -28,6 +28,7 @@ FixNVESphereKokkos<DeviceType>::FixNVESphereKokkos(LAMMPS *lmp, int narg, char *
   FixNVESphere(lmp, narg, arg)
 {
   kokkosable = 1;
+  fuse_integrate_flag = 1;
   atomKK = (AtomKokkos *)atom;
   execution_space = ExecutionSpaceFromDevice<DeviceType>::space;
 
@@ -164,6 +165,73 @@ void FixNVESphereKokkos<DeviceType>::final_integrate_item(const int i) const
   }
 }
 
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void FixNVESphereKokkos<DeviceType>::fused_integrate(int /*vflag*/)
+{
+  if (extra == DIPOLE)
+    atomKK->sync(execution_space, X_MASK | V_MASK | OMEGA_MASK| F_MASK | TORQUE_MASK | RMASS_MASK | RADIUS_MASK | MASK_MASK | MU_MASK);
+  else
+    atomKK->sync(execution_space, X_MASK | V_MASK | OMEGA_MASK| F_MASK | TORQUE_MASK | RMASS_MASK | RADIUS_MASK | MASK_MASK);
+
+  x = atomKK->k_x.view<DeviceType>();
+  v = atomKK->k_v.view<DeviceType>();
+  omega = atomKK->k_omega.view<DeviceType>();
+  f = atomKK->k_f.view<DeviceType>();
+  torque = atomKK->k_torque.view<DeviceType>();
+  mask = atomKK->k_mask.view<DeviceType>();
+  rmass = atomKK->k_rmass.view<DeviceType>();
+  radius = atomKK->k_radius.view<DeviceType>();
+  mu = atomKK->k_mu.view<DeviceType>();
+
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup) nlocal = atom->nfirst;
+
+  FixNVESphereKokkosFusedIntegrateFunctor<DeviceType> f(this);
+  Kokkos::parallel_for(nlocal,f);
+
+  if (extra == DIPOLE)
+    atomKK->modified(execution_space,  X_MASK | V_MASK | OMEGA_MASK | MU_MASK);
+  else
+    atomKK->modified(execution_space,  X_MASK | V_MASK | OMEGA_MASK);
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void FixNVESphereKokkos<DeviceType>::fused_integrate_item(const int i) const
+{
+  const double dtfrotate = dtf / inertia;
+
+  if (mask(i) & groupbit) {
+    const double dtfm = 2.0 * dtf / rmass(i);
+    v(i,0) += dtfm * f(i,0);
+    v(i,1) += dtfm * f(i,1);
+    v(i,2) += dtfm * f(i,2);
+    x(i,0) += dtv * v(i,0);
+    x(i,1) += dtv * v(i,1);
+    x(i,2) += dtv * v(i,2);
+
+    const double dtirotate = 2.0 * dtfrotate / (radius(i)*radius(i)*rmass(i));
+    omega(i,0) += dtirotate * torque(i,0);
+    omega(i,1) += dtirotate * torque(i,1);
+    omega(i,2) += dtirotate * torque(i,2);
+
+    if (extra == DIPOLE) {
+      const double g0 = mu(i,0) + dtv * (omega(i,1) * mu(i,2) - omega(i,2) * mu(i,1));
+      const double g1 = mu(i,1) + dtv * (omega(i,2) * mu(i,0) - omega(i,0) * mu(i,2));
+      const double g2 = mu(i,2) + dtv * (omega(i,0) * mu(i,1) - omega(i,1) * mu(i,0));
+      const double msq = g0*g0 + g1*g1 + g2*g2;
+      const double scale = mu(i,3)/sqrt(msq);
+      mu(i,0) = g0*scale;
+      mu(i,1) = g1*scale;
+      mu(i,2) = g2*scale;
+    }
+  }
+}
+
 namespace LAMMPS_NS {
 template class FixNVESphereKokkos<LMPDeviceType>;
 #ifdef LMP_KOKKOS_GPU

src/KOKKOS/fix_nve_sphere_kokkos.h

@@ -37,11 +37,14 @@ class FixNVESphereKokkos : public FixNVESphere {
     void init() override;
     void initial_integrate(int) override;
     void final_integrate() override;
+    void fused_integrate(int) override;
 
     KOKKOS_INLINE_FUNCTION
     void initial_integrate_item(const int i) const;
     KOKKOS_INLINE_FUNCTION
     void final_integrate_item(const int i) const;
+    KOKKOS_INLINE_FUNCTION
+    void fused_integrate_item(int) const;
 
   private:
     typename ArrayTypes<DeviceType>::t_x_array x;
@@ -77,6 +80,17 @@ struct FixNVESphereKokkosFinalIntegrateFunctor {
   }
 };
 
+template <class DeviceType>
+struct FixNVESphereKokkosFusedIntegrateFunctor {
+  typedef DeviceType device_type;
+  FixNVESphereKokkos<DeviceType> c;
+  FixNVESphereKokkosFusedIntegrateFunctor(FixNVESphereKokkos<DeviceType> *c_ptr): c(*c_ptr) { c.cleanup_copy(); }
+  KOKKOS_INLINE_FUNCTION
+  void operator()(const int i) const {
+    c.fused_integrate_item(i);
+  }
+};
+
 } // namespace LAMMPS_NS
 
 #endif // LMP_FIX_NVE_SPHERE_KOKKOS_H

src/KOKKOS/modify_kokkos.cpp

@@ -392,6 +392,24 @@ void ModifyKokkos::final_integrate()
   }
 }
 
+/* ----------------------------------------------------------------------
+   fused initial and final integrate call, only for relevant fixes
+------------------------------------------------------------------------- */
+
+void ModifyKokkos::fused_integrate(int vflag)
+{
+  for (int i = 0; i < n_final_integrate; i++) {
+    atomKK->sync(fix[list_final_integrate[i]]->execution_space,
+                 fix[list_final_integrate[i]]->datamask_read);
+    int prev_auto_sync = lmp->kokkos->auto_sync;
+    if (!fix[list_final_integrate[i]]->kokkosable) lmp->kokkos->auto_sync = 1;
+    fix[list_final_integrate[i]]->fused_integrate(vflag);
+    lmp->kokkos->auto_sync = prev_auto_sync;
+    atomKK->modified(fix[list_final_integrate[i]]->execution_space,
+                     fix[list_final_integrate[i]]->datamask_modify);
+  }
+}
+
 /* ----------------------------------------------------------------------
    end-of-timestep call, only for relevant fixes
    only call fix->end_of_step() on timesteps that are multiples of nevery
@@ -881,3 +899,22 @@ int ModifyKokkos::min_reset_ref()
   }
   return itmpall;
 }
+
+/* ----------------------------------------------------------------------
+   check if initial and final integrate can be fused
+------------------------------------------------------------------------- */
+
+int ModifyKokkos::check_fuse_integrate()
+{
+  int fuse_integrate_flag = 1;
+
+  for (int i = 0; i < n_initial_integrate; i++)
+    if (!fix[list_initial_integrate[i]]->fuse_integrate_flag)
+      fuse_integrate_flag = 0;
+
+  for (int i = 0; i < n_final_integrate; i++)
+    if (!fix[list_final_integrate[i]]->fuse_integrate_flag)
+      fuse_integrate_flag = 0;
+
+  return fuse_integrate_flag;
+}

src/KOKKOS/modify_kokkos.h

@@ -39,6 +39,7 @@ class ModifyKokkos : public Modify {
   void pre_reverse(int,int) override;
   void post_force(int) override;
   void final_integrate() override;
+  void fused_integrate(int) override;
   void end_of_step() override;
   double energy_couple() override;
   double energy_global() override;
@@ -69,6 +70,8 @@ class ModifyKokkos : public Modify {
   int min_dof() override;
   int min_reset_ref() override;
 
+  int check_fuse_integrate();
+
  protected:
 
 };

src/KOKKOS/pair_kokkos.h

@@ -280,6 +280,12 @@ struct PairComputeFunctor  {
       const X_FLOAT ztmp = c.x(i,2);
       const int itype = c.type(i);
 
+      Kokkos::single(Kokkos::PerThread(team), [&] (){
+        f(i,0) = 0.0;
+        f(i,1) = 0.0;
+        f(i,2) = 0.0;
+      });
+
       const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i);
       const int jnum = list.d_numneigh[i];
 
@@ -337,6 +343,12 @@ struct PairComputeFunctor  {
       const int itype = c.type(i);
       const F_FLOAT qtmp = c.q(i);
 
+      Kokkos::single(Kokkos::PerThread(team), [&] (){
+        f(i,0) = 0.0;
+        f(i,1) = 0.0;
+        f(i,2) = 0.0;
+      });
+
       const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i);
       const int jnum = list.d_numneigh[i];
 
@@ -399,6 +411,12 @@ struct PairComputeFunctor  {
       const X_FLOAT ztmp = c.x(i,2);
       const int itype = c.type(i);
 
+      Kokkos::single(Kokkos::PerThread(team), [&] (){
+        f(i,0) = 0.0;
+        f(i,1) = 0.0;
+        f(i,2) = 0.0;
+      });
+
       const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i);
       const int jnum = list.d_numneigh[i];
 
@@ -495,6 +513,12 @@ struct PairComputeFunctor  {
       const int itype = c.type(i);
       const F_FLOAT qtmp = c.q(i);
 
+      Kokkos::single(Kokkos::PerThread(team), [&] (){
+        f(i,0) = 0.0;
+        f(i,1) = 0.0;
+        f(i,2) = 0.0;
+      });
+
       const AtomNeighborsConst neighbors_i = list.get_neighbors_const(i);
       const int jnum = list.d_numneigh[i];
 
@@ -743,6 +767,8 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename std::enable_if<(NEIG
       fpair->lmp->kokkos->neigh_thread = 1;
 
   if (fpair->lmp->kokkos->neigh_thread) {
+    fpair->fuse_force_clear_flag = 1;
+
     int vector_length = 8;
     int atoms_per_team = 32;
 

src/KOKKOS/verlet_kokkos.cpp

@@ -27,7 +27,7 @@
 #include "kspace.h"
 #include "output.h"
 #include "update.h"
-#include "modify.h"
+#include "modify_kokkos.h"
 #include "timer.h"
 #include "memory_kokkos.h"
 #include "kokkos.h"
@@ -276,6 +276,9 @@ void VerletKokkos::run(int n)
 
   lmp->kokkos->auto_sync = 0;
 
+  fuse_integrate = 0;
+  fuse_force_clear = 0;
+
   if (atomKK->sortfreq > 0) sortflag = 1;
   else sortflag = 0;
 
@@ -296,7 +299,8 @@ void VerletKokkos::run(int n)
     // initial time integration
 
     timer->stamp();
-    modify->initial_integrate(vflag);
+    if (!fuse_integrate)
+      modify->initial_integrate(vflag);
     if (n_post_integrate) modify->post_integrate();
     timer->stamp(Timer::MODIFY);
 
@@ -357,12 +361,17 @@ void VerletKokkos::run(int n)
       }
     }
 
+    // check if kernels can be fused, must come after initial_integrate
+
+    fuse_check(i,n);
+
     // force computations
     // important for pair to come before bonded contributions
     // since some bonded potentials tally pairwise energy/virial
     // and Pair:ev_tally() needs to be called before any tallying
 
-    force_clear();
+    if (!fuse_force_clear)
+      force_clear();
 
     timer->stamp();
 
@@ -494,7 +503,10 @@ void VerletKokkos::run(int n)
     // force modifications, final time integration, diagnostics
 
     if (n_post_force) modify->post_force(vflag);
-    modify->final_integrate();
+
+    if (fuse_integrate) modify->fused_integrate(vflag);
+    else modify->final_integrate();
+
     if (n_end_of_step) modify->end_of_step();
     timer->stamp(Timer::MODIFY);
 
@@ -593,3 +605,35 @@ void VerletKokkos::force_clear()
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   check if can fuse force_clear() with pair compute()
+   Requirements:
+   - no pre_force fixes
+   - no torques, SPIN forces, or includegroup set
+   - pair compute() must be called
+   - pair_style must support fusing
+
+   check if can fuse initial_integrate() with final_integrate()
+   Requirements:
+   - no end_of_step fixes
+   - not on last or output step
+   - no timers to break out of loop
+   - integrate fix style must support fusing
+------------------------------------------------------------------------- */
+
+void VerletKokkos::fuse_check(int i, int n)
+{
+  fuse_force_clear = 1;
+  if (modify->n_pre_force) fuse_force_clear = 0;
+  else if (torqueflag || extraflag || neighbor->includegroup) fuse_force_clear = 0;
+  else if (!force->pair || !pair_compute_flag) fuse_force_clear = 0;
+  else if (!force->pair->fuse_force_clear_flag) fuse_force_clear = 0;
+
+  fuse_integrate = 1;
+  if (modify->n_end_of_step) fuse_integrate = 0;
+  else if (i == n-1) fuse_integrate = 0;
+  else if (update->ntimestep == output->next) fuse_integrate = 0;
+  else if (timer->has_timeout()) fuse_integrate = 0;
+  else if (!((ModifyKokkos*)modify)->check_fuse_integrate()) fuse_integrate = 0;
+}

src/KOKKOS/verlet_kokkos.h

@@ -46,6 +46,9 @@ class VerletKokkos : public Verlet {
 
  protected:
   DAT::t_f_array f_merge_copy,f;
+  int fuse_force_clear,fuse_integrate;
+
+  void fuse_check(int, int);
 };
 }
 

src/compute_count_type.cpp

@@ -0,0 +1,400 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "compute_count_type.h"
+
+#include "atom.h"
+#include "domain.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+
+enum { ATOM, BOND, ANGLE, DIHEDRAL, IMPROPER };
+
+/* ---------------------------------------------------------------------- */
+
+ComputeCountType::ComputeCountType(LAMMPS *lmp, int narg, char **arg) :
+    Compute(lmp, narg, arg), count(nullptr), bcount_me(nullptr), bcount(nullptr)
+{
+  if (narg != 4) error->all(FLERR, "Incorrect number of args for compute count/type command");
+
+  // process args
+
+  if (strcmp(arg[3], "atom") == 0)
+    mode = ATOM;
+  else if (strcmp(arg[3], "bond") == 0)
+    mode = BOND;
+  else if (strcmp(arg[3], "angle") == 0)
+    mode = ANGLE;
+  else if (strcmp(arg[3], "dihedral") == 0)
+    mode = DIHEDRAL;
+  else if (strcmp(arg[3], "improper") == 0)
+    mode = IMPROPER;
+  else
+    error->all(FLERR, "Invalid compute count/type keyword {}", arg[3]);
+
+  // error check
+
+  if (mode == BOND && !atom->nbondtypes)
+    error->all(FLERR, "Compute count/type bond command with no bonds defined");
+  if (mode == ANGLE && !atom->nangletypes)
+    error->all(FLERR, "Compute count/type bond command with no angles defined");
+  if (mode == DIHEDRAL && !atom->ndihedraltypes)
+    error->all(FLERR, "Compute count/type dihedral command with no dihedrals defined");
+  if (mode == IMPROPER && !atom->nimpropertypes)
+    error->all(FLERR, "Compute count/type improper command with no impropers defined");
+
+  // set vector lengths
+
+  if (mode == ATOM) {
+    vector_flag = 1;
+    size_vector = atom->ntypes;
+    extvector = 1;
+  } else if (mode == BOND) {
+    scalar_flag = vector_flag = 1;
+    size_vector = atom->nbondtypes;
+    extscalar = 1;
+    extvector = 1;
+  } else if (mode == ANGLE) {
+    vector_flag = 1;
+    size_vector = atom->nangletypes;
+    extvector = 1;
+  } else if (mode == DIHEDRAL) {
+    vector_flag = 1;
+    size_vector = atom->ndihedraltypes;
+    extvector = 1;
+  } else if (mode == IMPROPER) {
+    vector_flag = 1;
+    size_vector = atom->nimpropertypes;
+    extvector = 1;
+  }
+
+  // output vector
+
+  vector = new double[size_vector];
+
+  // work vectors
+
+  count = new int[size_vector];
+  bcount_me = new bigint[size_vector];
+  bcount = new bigint[size_vector];
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeCountType::~ComputeCountType()
+{
+  delete[] vector;
+
+  delete[] count;
+  delete[] bcount_me;
+  delete[] bcount;
+}
+
+/* ----------------------------------------------------------------------
+   only invoked for mode = BOND to count broken bonds
+   broken bonds have bond_type = 0
+---------------------------------------------------------------------- */
+
+double ComputeCountType::compute_scalar()
+{
+  invoked_scalar = update->ntimestep;
+
+  int *num_bond = atom->num_bond;
+  int **bond_type = atom->bond_type;
+  int nlocal = atom->nlocal;
+
+  // count broken bonds with bond_type = 0
+  // ignore group setting since 2 atoms in a broken bond
+  //   can be arbitrarily far apart
+
+  int count = 0;
+  for (int i = 0; i < nlocal; i++) {
+    int nbond = num_bond[i];
+    for (int m = 0; m < nbond; m++)
+      if (bond_type[i][m] == 0) count++;
+  }
+
+  // sum across procs as bigint, then convert to double
+  // correct for double counting if newton_bond off
+
+  bigint bcount = 0;
+  bigint bcount_me = count;
+  MPI_Allreduce(&bcount_me, &bcount, 1, MPI_LMP_BIGINT, MPI_SUM, world);
+  if (force->newton_bond == 0) bcount /= 2;
+
+  if (bcount > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow");
+  scalar = bcount;
+  return scalar;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeCountType::compute_vector()
+{
+  invoked_vector = update->ntimestep;
+
+  int nvec;
+
+  if (mode == ATOM)
+    nvec = count_atoms();
+  else if (mode == BOND)
+    nvec = count_bonds();
+  else if (mode == ANGLE)
+    nvec = count_angles();
+  else if (mode == DIHEDRAL)
+    nvec = count_dihedrals();
+  else if (mode == IMPROPER)
+    nvec = count_impropers();
+
+  // sum across procs as bigint, then convert to double
+  // correct for multiple counting if newton_bond off
+
+  for (int m = 0; m < nvec; m++) bcount_me[m] = count[m];
+  MPI_Allreduce(bcount_me, bcount, nvec, MPI_LMP_BIGINT, MPI_SUM, world);
+
+  if (force->newton_bond == 0) {
+    if (mode == BOND)
+      for (int m = 0; m < nvec; m++) bcount[m] /= 2;
+    else if (mode == ANGLE)
+      for (int m = 0; m < nvec; m++) bcount[m] /= 3;
+    if (mode == DIHEDRAL || mode == IMPROPER)
+      for (int m = 0; m < nvec; m++) bcount[m] /= 4;
+  }
+
+  for (int m = 0; m < nvec; m++)
+    if (bcount[m] > MAXDOUBLEINT) error->all(FLERR, "Compute count/type overflow");
+  for (int m = 0; m < nvec; m++) vector[m] = bcount[m];
+}
+
+/* ----------------------------------------------------------------------
+   count atoms by type
+   atom must be in group to be counted
+---------------------------------------------------------------------- */
+
+int ComputeCountType::count_atoms()
+{
+  int *type = atom->type;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int ntypes = atom->ntypes;
+
+  for (int m = 0; m < ntypes; m++) count[m] = 0;
+  for (int i = 0; i < nlocal; i++)
+    if (mask[i] & groupbit) count[type[i] - 1]++;
+
+  return ntypes;
+}
+
+/* ----------------------------------------------------------------------
+   count bonds by type
+   both atoms in bond must be in group to be counted
+   skip type = 0 bonds, they are counted by compute_scalar()
+   bond types can be negative, count them as if positive
+---------------------------------------------------------------------- */
+
+int ComputeCountType::count_bonds()
+{
+  tagint **bond_atom = atom->bond_atom;
+  int **bond_type = atom->bond_type;
+  int *num_bond = atom->num_bond;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int nbondtypes = atom->nbondtypes;
+
+  int flag = 0;
+  for (int m = 0; m < nbondtypes; m++) count[m] = 0;
+
+  for (int i = 0; i < nlocal; i++) {
+    int nbond = num_bond[i];
+    for (int m = 0; m < nbond; m++) {
+      int itype = bond_type[i][m];
+      if (itype == 0) continue;
+
+      int j = atom->map(bond_atom[i][m]);
+      if (j < 0) {
+        flag = 1;
+        continue;
+      }
+
+      if ((mask[i] & groupbit) && (mask[j] & groupbit)) {
+        if (itype > 0)
+          count[itype - 1]++;
+        else
+          count[-itype - 1]++;
+      }
+    }
+  }
+
+  int flagany;
+  MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world);
+  if (flagany) error->all(FLERR, "Missing bond atom in compute count/type");
+
+  return nbondtypes;
+}
+
+/* ----------------------------------------------------------------------
+   count angles by type
+   all 3 atoms in angle must be in group to be counted
+   angle types can be negative, count them as if positive
+---------------------------------------------------------------------- */
+
+int ComputeCountType::count_angles()
+{
+  tagint **angle_atom1 = atom->angle_atom1;
+  tagint **angle_atom2 = atom->angle_atom2;
+  tagint **angle_atom3 = atom->angle_atom3;
+  int **angle_type = atom->angle_type;
+  int *num_angle = atom->num_angle;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int nangletypes = atom->nangletypes;
+
+  int flag = 0;
+  for (int m = 0; m < nangletypes; m++) count[m] = 0;
+
+  for (int i = 0; i < nlocal; i++) {
+    int nangle = num_angle[i];
+    for (int m = 0; m < nangle; m++) {
+      int itype = angle_type[i][m];
+
+      int j1 = atom->map(angle_atom1[i][m]);
+      int j2 = atom->map(angle_atom2[i][m]);
+      int j3 = atom->map(angle_atom3[i][m]);
+      if (j1 < 0 || j2 < 0 || j3 < 0) {
+        flag = 1;
+        continue;
+      }
+
+      if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit)) {
+        if (itype > 0)
+          count[itype - 1]++;
+        else if (itype < 0)
+          count[-itype - 1]++;
+      }
+    }
+  }
+
+  int flagany;
+  MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world);
+  if (flagany) error->all(FLERR, "Missing angle atom in compute count/type");
+
+  return nangletypes;
+}
+
+/* ----------------------------------------------------------------------
+   count dihedrals by type
+   all 4 atoms in dihedral must be in group to be counted
+   dihedral types can be negative, count them as if positive
+---------------------------------------------------------------------- */
+
+int ComputeCountType::count_dihedrals()
+{
+  tagint **dihedral_atom1 = atom->dihedral_atom1;
+  tagint **dihedral_atom2 = atom->dihedral_atom2;
+  tagint **dihedral_atom3 = atom->dihedral_atom3;
+  tagint **dihedral_atom4 = atom->dihedral_atom4;
+  int **dihedral_type = atom->dihedral_type;
+  int *num_dihedral = atom->num_dihedral;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int ndihedraltypes = atom->ndihedraltypes;
+
+  int flag = 0;
+  for (int m = 0; m < ndihedraltypes; m++) count[m] = 0;
+
+  for (int i = 0; i < nlocal; i++) {
+    int ndihedral = num_dihedral[i];
+    for (int m = 0; m < ndihedral; m++) {
+      int itype = dihedral_type[i][m];
+
+      int j1 = atom->map(dihedral_atom1[i][m]);
+      int j2 = atom->map(dihedral_atom2[i][m]);
+      int j3 = atom->map(dihedral_atom3[i][m]);
+      int j4 = atom->map(dihedral_atom4[i][m]);
+      if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) {
+        flag = 1;
+        continue;
+      }
+
+      if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) &&
+          (mask[j4] & groupbit)) {
+        if (itype > 0)
+          count[itype - 1]++;
+        else if (itype < 0)
+          count[-itype - 1]++;
+      }
+    }
+  }
+
+  int flagany;
+  MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world);
+  if (flagany) error->all(FLERR, "Missing dihedral atom in compute count/type");
+
+  return ndihedraltypes;
+}
+
+/* ----------------------------------------------------------------------
+   count impropers by type
+   all 4 atoms in improper must be in group to be counted
+   improper types can be negative, count them as if positive
+---------------------------------------------------------------------- */
+
+int ComputeCountType::count_impropers()
+{
+  tagint **improper_atom1 = atom->improper_atom1;
+  tagint **improper_atom2 = atom->improper_atom2;
+  tagint **improper_atom3 = atom->improper_atom3;
+  tagint **improper_atom4 = atom->improper_atom4;
+  int **improper_type = atom->improper_type;
+  int *num_improper = atom->num_improper;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int nimpropertypes = atom->nimpropertypes;
+
+  int flag = 0;
+  for (int m = 0; m < nimpropertypes; m++) count[m] = 0;
+
+  for (int i = 0; i < nlocal; i++) {
+    int nimproper = num_improper[i];
+    for (int m = 0; m < nimproper; m++) {
+      int itype = improper_type[i][m];
+
+      int j1 = atom->map(improper_atom1[i][m]);
+      int j2 = atom->map(improper_atom2[i][m]);
+      int j3 = atom->map(improper_atom3[i][m]);
+      int j4 = atom->map(improper_atom4[i][m]);
+      if (j1 < 0 || j2 < 0 || j3 < 0 || j4 < 0) {
+        flag = 1;
+        continue;
+      }
+
+      if ((mask[j1] & groupbit) && (mask[j2] & groupbit) && (mask[j3] & groupbit) &&
+          (mask[j4] & groupbit)) {
+        if (itype > 0)
+          count[itype - 1]++;
+        else if (itype < 0)
+          count[-itype - 1]++;
+      }
+    }
+  }
+
+  int flagany;
+  MPI_Allreduce(&flag, &flagany, 1, MPI_INT, MPI_SUM, world);
+  if (flagany) error->all(FLERR, "Missing improper atom in compute count/type");
+
+  return nimpropertypes;
+}

src/compute_count_type.h

@@ -0,0 +1,52 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   https://www.lammps.org/, Sandia National Laboratories
+   LAMMPS development team: developers@lammps.org
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef COMPUTE_CLASS
+// clang-format off
+ComputeStyle(count/type,ComputeCountType);
+// clang-format on
+#else
+
+#ifndef LMP_COMPUTE_COMPUTE_TYPE_H
+#define LMP_COMPUTE_COMPUTE_TYPE_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeCountType : public Compute {
+ public:
+  ComputeCountType(class LAMMPS *, int, char **);
+  ~ComputeCountType() override;
+  void init() override {}
+  double compute_scalar() override;
+  void compute_vector() override;
+
+ protected:
+  int mode;
+
+  int *count;
+  bigint *bcount_me;
+  bigint *bcount;
+
+  int count_atoms();
+  int count_bonds();
+  int count_angles();
+  int count_dihedrals();
+  int count_impropers();
+};
+
+}    // namespace LAMMPS_NS
+
+#endif
+#endif

src/fix.cpp

@@ -102,15 +102,15 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) :
   vflag_atom = cvflag_atom = 0;
   centroidstressflag = CENTROID_SAME;
 
-  // KOKKOS per-fix data masks
+  // KOKKOS package
 
   execution_space = Host;
   datamask_read = ALL_MASK;
   datamask_modify = ALL_MASK;
 
-  kokkosable = 0;
+  kokkosable = copymode = 0;
   forward_comm_device = exchange_comm_device = sort_device = 0;
-  copymode = 0;
+  fuse_integrate_flag = 0;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix.h

@@ -127,11 +127,12 @@ class Fix : protected Pointers {
 
   int restart_reset;    // 1 if restart just re-initialized fix
 
-  // KOKKOS host/device flag and data masks
+  // KOKKOS flags and variables
 
   int kokkosable;             // 1 if Kokkos fix
   int forward_comm_device;    // 1 if forward comm on Device
   int exchange_comm_device;   // 1 if exchange comm on Device
+  int fuse_integrate_flag;    // 1 if can fuse initial integrate with final integrate
   int sort_device;            // 1 if sort on Device
   ExecutionSpace execution_space;
   unsigned int datamask_read, datamask_modify;
@@ -161,6 +162,7 @@ class Fix : protected Pointers {
   virtual void pre_reverse(int, int) {}
   virtual void post_force(int) {}
   virtual void final_integrate() {}
+  virtual void fused_integrate(int) {}
   virtual void end_of_step() {}
   virtual void post_run() {}
   virtual void write_restart(FILE *) {}

src/lmptype.h

@@ -95,6 +95,7 @@ typedef int64_t bigint;
 #define MAXSMALLINT INT_MAX
 #define MAXTAGINT INT_MAX
 #define MAXBIGINT INT64_MAX
+#define MAXDOUBLEINT 9007199254740992   // 2^53
 
 #define MPI_LMP_TAGINT MPI_INT
 #define MPI_LMP_IMAGEINT MPI_INT
@@ -132,6 +133,7 @@ typedef int64_t bigint;
 #define MAXSMALLINT INT_MAX
 #define MAXTAGINT INT64_MAX
 #define MAXBIGINT INT64_MAX
+#define MAXDOUBLEINT 9007199254740992   // 2^53
 
 #define MPI_LMP_TAGINT MPI_LL
 #define MPI_LMP_IMAGEINT MPI_LL
@@ -168,6 +170,7 @@ typedef int bigint;
 #define MAXSMALLINT INT_MAX
 #define MAXTAGINT INT_MAX
 #define MAXBIGINT INT_MAX
+#define MAXDOUBLEINT INT_MAX
 
 #define MPI_LMP_TAGINT MPI_INT
 #define MPI_LMP_IMAGEINT MPI_INT

src/modify.h

@@ -69,6 +69,7 @@ class Modify : protected Pointers {
   virtual void pre_reverse(int, int);
   virtual void post_force(int);
   virtual void final_integrate();
+  virtual void fused_integrate(int) {}
   virtual void end_of_step();
   virtual double energy_couple();
   virtual double energy_global();

src/pair.cpp

@@ -116,15 +116,14 @@ Pair::Pair(LAMMPS *lmp) :
   nondefault_history_transfer = 0;
   beyond_contact = 0;
 
-  // KOKKOS per-fix data masks
+  // KOKKOS package
 
   execution_space = Host;
   datamask_read = ALL_MASK;
   datamask_modify = ALL_MASK;
 
-  kokkosable = 0;
+  kokkosable = copymode = 0;
-  reverse_comm_device = 0;
+  reverse_comm_device = fuse_force_clear_flag = 0;
-  copymode = 0;
 }
 
 /* ---------------------------------------------------------------------- */

src/pair.h

@@ -119,12 +119,13 @@ class Pair : protected Pointers {
 
   int beyond_contact, nondefault_history_transfer;    // for granular styles
 
-  // KOKKOS host/device flag and data masks
+  // KOKKOS flags and variables
 
   ExecutionSpace execution_space;
   unsigned int datamask_read, datamask_modify;
   int kokkosable;             // 1 if Kokkos pair
   int reverse_comm_device;    // 1 if reverse comm on Device
+  int fuse_force_clear_flag;   // 1 if can fuse force clear with force compute
 
   Pair(class LAMMPS *);
   ~Pair() override;

src/timer.h

@@ -63,6 +63,7 @@ class Timer : protected Pointers {
   bool has_normal() const { return (_level >= NORMAL); }
   bool has_full() const { return (_level >= FULL); }
   bool has_sync() const { return (_sync != OFF); }
+  bool has_timeout() const { return (_timeout >= 0.0); }
 
   // flag if wallclock time is expired
   bool is_timeout() const { return (_timeout == 0.0); }

tools/phonon/CMakeLists.spglib

@@ -0,0 +1,54 @@
+cmake_minimum_required(VERSION 3.10)
+project(spglib C)
+set(CMAKE_MACOSX_RPATH 1)
+set(CMAKE_C_FLAGS_RELEASE "-Wall -O2")
+set(CMAKE_C_FLAGS_DEBUG "-g -DSPGDEBUG -DSPGWARNING")
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+endif(NOT CMAKE_BUILD_TYPE)
+message(STATUS "Build type: ${CMAKE_BUILD_TYPE}")
+
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+
+# Version numbers
+file(READ ${PROJECT_SOURCE_DIR}/src/version.h version_file)
+string(REGEX MATCH "SPGLIB_MAJOR_VERSION ([0-9]+)" spglib_major_version ${version_file})
+set(spglib_major_version ${CMAKE_MATCH_1})
+string(REGEX MATCH "SPGLIB_MINOR_VERSION ([0-9]+)" spglib_minor_version ${version_file})
+set(spglib_minor_version ${CMAKE_MATCH_1})
+string(REGEX MATCH "SPGLIB_MICRO_VERSION ([0-9]+)" spglib_micro_version ${version_file})
+set(spglib_micro_version ${CMAKE_MATCH_1})
+set(serial "${spglib_major_version}.${spglib_minor_version}.${spglib_micro_version}")
+set(soserial "1")
+
+# Source code
+include_directories("${PROJECT_SOURCE_DIR}/src")
+set(SOURCES ${PROJECT_SOURCE_DIR}/src/arithmetic.c
+            ${PROJECT_SOURCE_DIR}/src/cell.c
+            ${PROJECT_SOURCE_DIR}/src/debug.c
+            ${PROJECT_SOURCE_DIR}/src/delaunay.c
+            ${PROJECT_SOURCE_DIR}/src/determination.c
+            ${PROJECT_SOURCE_DIR}/src/hall_symbol.c
+            ${PROJECT_SOURCE_DIR}/src/kgrid.c
+            ${PROJECT_SOURCE_DIR}/src/kpoint.c
+            ${PROJECT_SOURCE_DIR}/src/mathfunc.c
+            ${PROJECT_SOURCE_DIR}/src/niggli.c
+            ${PROJECT_SOURCE_DIR}/src/overlap.c
+            ${PROJECT_SOURCE_DIR}/src/pointgroup.c
+            ${PROJECT_SOURCE_DIR}/src/primitive.c
+            ${PROJECT_SOURCE_DIR}/src/refinement.c
+            ${PROJECT_SOURCE_DIR}/src/site_symmetry.c
+            ${PROJECT_SOURCE_DIR}/src/sitesym_database.c
+            ${PROJECT_SOURCE_DIR}/src/spacegroup.c
+            ${PROJECT_SOURCE_DIR}/src/spg_database.c
+            ${PROJECT_SOURCE_DIR}/src/spglib.c
+            ${PROJECT_SOURCE_DIR}/src/spin.c
+            ${PROJECT_SOURCE_DIR}/src/symmetry.c)
+
+add_library(symspg STATIC ${SOURCES})
+install(TARGETS symspg ARCHIVE DESTINATION ${CMAKE_INSTALL_PREFIX}/lib)
+
+# Header file
+install(FILES ${PROJECT_SOURCE_DIR}/src/spglib.h DESTINATION ${CMAKE_INSTALL_PREFIX}/include)
+

tools/phonon/CMakeLists.txt

@@ -0,0 +1,120 @@
+
+# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL),
+# macOS, MSVC 2019 (=Version 16)
+cmake_minimum_required(VERSION 3.10)
+
+# set timestamp of downloaded files to that of archive
+if(POLICY CMP0135)
+  cmake_policy(SET CMP0135 NEW)
+endif()
+
+# set up project
+set(PHANA_MINOR_VERSION 48)
+project(phonon VERSION ${PHANA_MINOR_VERSION}
+  DESCRIPTION "Fix phonon post-processor"
+  LANGUAGES CXX C)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE RelWithDebInfo)
+endif()
+
+# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions,
+# padding and Spectre mitigation
+if(MSVC)
+  add_compile_options(/wd4244)
+  add_compile_options(/wd4267)
+  add_compile_options(/wd4711)
+  add_compile_options(/wd4820)
+  add_compile_options(/wd5045)
+  add_compile_definitions(_CRT_SECURE_NO_WARNINGS)
+endif()
+
+set(CMAKE_PROJECT_VERSION ${PHANA_MINOR_VERSION})
+configure_file(version.h.in version.h @ONLY)
+add_executable(phana
+  main.cpp
+  disp.cpp
+  dynmat.cpp
+  green.cpp
+  input.cpp
+  interpolate.cpp
+  kpath.cpp
+  memory.cpp
+  phonon.cpp
+  phonopy.cpp
+  qnodes.cpp
+  timer.cpp
+)
+target_include_directories(phana PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>)
+
+if(NOT LAMMPS_DIR)
+  set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/../../cmake/Modules)
+  set(LAMMPS_THIRDPARTY_URL "https://download.lammps.org/thirdparty")
+endif()
+find_package(FFTW3)
+if(FFTW3_FOUND)
+  target_compile_definitions(phana PRIVATE FFTW3)
+  target_link_libraries(phana PRIVATE FFTW3::FFTW3)
+endif()
+
+# build bundeled libraries
+add_subdirectory(tricubic)
+
+# standalone build must build our own version of linalg
+if(NOT LAMMPS_DIR)
+  if(NOT USE_INTERNAL_LINALG)
+    find_package(LAPACK)
+    find_package(BLAS)
+  endif()
+  if(NOT LAPACK_FOUND OR NOT BLAS_FOUND OR USE_INTERNAL_LINALG)
+    file(GLOB LINALG_SOURCES CONFIGURE_DEPENDS ${CMAKE_CURRENT_SOURCE_DIR}/../../lib/linalg/[^.]*.cpp)
+    add_library(linalg STATIC ${LINALG_SOURCES})
+    set(BLAS_LIBRARIES "$<TARGET_FILE:linalg>")
+    set(LAPACK_LIBRARIES "$<TARGET_FILE:linalg>")
+  else()
+    list(APPEND LAPACK_LIBRARIES ${BLAS_LIBRARIES})
+  endif()
+endif()
+
+option(USE_SPGLIB "Download and use spglib for phonon DOS and other optional properties" ON)
+if(USE_SPGLIB)
+  set(SPGLIB_URL "https://github.com/spglib/spglib/archive/refs/tags/v1.11.2.1.tar.gz" CACHE STRING "URL for spglib v1.x tarball")
+  set(SPGLIB_MD5 "3089782bc85b5034dd4765a18ee70bc7" CACHE STRING "MD5 checksum for spglib tarball")
+  mark_as_advanced(SPGLIB_URL)
+  mark_as_advanced(SPGLIB_MD5)
+  include(LAMMPSUtils)
+  GetFallbackURL(SPGLIB_URL SPGLIB_FALLBACK)
+
+  include(ExternalProject)
+  ExternalProject_Add(spglib_build
+    URL ${SPGLIB_URL} ${SPGLIB_FALLBACK}
+    URL_MD5 ${SPGLIB_MD5}
+    PREFIX ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext
+    CMAKE_ARGS -DCMAKE_INSTALL_PREFIX=${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext
+    -DCMAKE_C_COMPILER=${CMAKE_C_COMPILER}
+    -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
+    -DCMAKE_MAKE_PROGRAM=${CMAKE_MAKE_PROGRAM}
+    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE}
+    -DCMAKE_POSITION_INDEPENDENT_CODE=ON
+    UPDATE_COMMAND ${CMAKE_COMMAND} -E copy_if_different ${CMAKE_CURRENT_SOURCE_DIR}/CMakeLists.spglib ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build/CMakeLists.txt
+    INSTALL_COMMAND ${CMAKE_COMMAND} --build ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/src/spglib_build-build --target install
+    BUILD_BYPRODUCTS "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}"
+  )
+
+  # workaround for older CMake versions
+  file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib)
+  file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include)
+
+  add_library(SPGLIB::SYMSPG UNKNOWN IMPORTED)
+  add_dependencies(SPGLIB::SYMSPG spglib_build)
+  set_target_properties(SPGLIB::SYMSPG PROPERTIES
+    IMPORTED_LOCATION "${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/lib/${CMAKE_STATIC_LIBRARY_PREFIX}symspg${CMAKE_STATIC_LIBRARY_SUFFIX}"
+    INTERFACE_INCLUDE_DIRECTORIES ${CMAKE_CURRENT_BINARY_DIR}/spglib_build_ext/include
+  )
+    
+  target_compile_definitions(phana PRIVATE UseSPG)
+  target_link_libraries(phana PRIVATE SPGLIB::SYMSPG)
+endif()
+
+target_link_libraries(phana PRIVATE tricubic ${LAPACK_LIBRARIES})
+install(TARGETS phana EXPORT LAMMPS_Targets DESTINATION ${CMAKE_INSTALL_BINDIR})

tools/phonon/Makefile

@@ -1,67 +0,0 @@
-.SUFFIXES : .o .cpp
-# compiler and flags
-CC     = g++ -Wno-unused-result
-LINK   = $(CC) -static
-CFLAGS = -O3 $(DEBUG) $(UFLAG)
-
-OFLAGS = -O3 $(DEBUG)
-INC    = $(LPKINC) $(TCINC) $(SPGINC) $(FFTINC)
-LIB    = $(LPKLIB) $(TCLIB) $(SPGLIB) $(FFTLIB)
-
-# cLapack library needed
-LPKINC = -I/opt/clapack/3.2.1/include
-LPKLIB = -L/opt/clapack/3.2.1/lib -lclapack -lblas -lf2c -lm
-
-# Tricubic library needed
-TCINC = -I/opt/tricubic/1.0/include
-TCLIB = -L/opt/tricubic/1.0/lib -ltricubic
-
-# spglib, used to get the irreducible q-points
-# if SFLAG is not set, spglib won't be used.
-SFLAG  = -DUseSPG
-SPGINC = -I/opt/spglib/1.9.7/include/spglib
-SPGLIB = -L/opt/spglib/1.9.7/lib -lsymspg
-
-# FFTW 3， used to deduce the force constants in real space
-# if FFLAG is not set, fftw won't be used.
-FFLAG  = -DFFTW3
-FFTINC = -I/opt/fftw/3.3.7/include
-FFTLIB = -L/opt/fftw/3.3.7/lib -lfftw3
-
-# Debug flags
-# DEBUG = -g -DDEBUG
-UFLAG = $(SFLAG) $(FFLAG)
-
-#====================================================================
-ROOT   = phana
-# executable name
-EXE    = $(ROOT)
-#====================================================================
-# source and rules
-SRC = $(wildcard *.cpp)
-OBJ = $(SRC:.cpp=.o)
-
-#====================================================================
-all:  ${EXE}
-
-${EXE}: $(OBJ)
-	$(LINK) $(OFLAGS) $(OBJ) $(LIB) -o $@
-
-clean: 
-	rm -f *.o *~ *.mod ${EXE}
-
-tar:
-	rm -f ${ROOT}.tar.gz; tar -czvf ${ROOT}.tar.gz *.cpp  *.h Makefile README
-
-ver:
-	@echo "#define VERSION `git log|grep '^commit'|wc -l`" > version.h
-
-#====================================================================
-.f.o:
-	$(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $<
-.f90.o:
-	$(FC) $(FFLAGS) $(FREE) $(MPI) ${INC} -c $<
-.c.o:
-	$(CC) $(CFLAGS) -c $<
-.cpp.o:
-	$(CC) $(CFLAGS) $(INC) -c $<

tools/phonon/README

@@ -5,34 +5,38 @@
    analyse the phonon related information.
 #-------------------------------------------------------------------------------
 1. Dependencies
-   The clapack library is needed to solve the eigen problems,
+   The ZHEEVD LAPACK function is needed to solve the eigen problems.
-   which could be downloaded from:
+   A C++ compilable version based on CLAPACK is included in the linalg folder
-   http://www.netlib.org/clapack/
+   and will be automatically built.
-   
+
-   The tricubic library is also needed to do tricubic interpolations,
+   The tricubic library is also needed to do tricubic interpolations.
-   which could now be obtained from:
+   A copy is included and will be automatically built.
-   https://github.com/nbigaouette/libtricubic/
+
-   
    The spglib is optionally needed, enabling one to evaluate the
    phonon density of states or vibrational thermal properties
    using only the irreducible q-points in the first Brillouin zone,
    as well as to evaluate the phonon dispersion curvers with the
-   automatic mode. Currently, the 1.8.3 version of spglib is used.
+   automatic mode. Currently, version 1.11.2.1 of spglib is used.
-   It can be obtained from:
+   It is automatically downloaded and compiled unless the -DUSE_SPGLIB=off
-   http://spglib.sourceforge.net/
+   variable is set during CMake configuration.
 
    FFTW 3 might also be needed if you would like to interface with
    phonopy: necessary input files for phonopy will be prepared so
    that you can make use of the functions provided by phonopy.
-   FFTW 3 can be downloaded from:
+   It is autodetected and used if available.
-   http://www.fftw.org
+
-   
+   FFTW 3 can be downloaded from: http://www.fftw.org
+
 2. Compilation
-   To compile the code, one needs therefore to install the above
+   To compile the code, one needs to have CMake version 3.16
-   libraries and set the paths correctly in the Makefile.
+   or later installed.
-   Once this is done, by typing
+
-     make
+   The CMake configuration is done with:
-   will yield the executable "phana".
+     cmake -S . -B build
+   And compilation then performed with:
+     cmake --build build
+   The phana (or phana.exe) executable is then available in
+   the "build" folder
 
 3. Unit system
    The units of the output frequencies by this code is THz for
@@ -46,6 +50,18 @@
 5. Bug report
    If any bug found, please drop a line to: konglt(at)sjtu.edu.cn
 
+6. Precompiled executable
+   The "precompiled" folder contains a precompiled and statically
+   linked Linux executable for x86_64 CPUs. It should work on *any*
+   Linux machine with using the x86_64 architecture. It includes
+   spglib support but not fftw3.
+
+7. Portability
+   Build and use of phana has been successfully tested on:
+    - Fedora Linux 38 using GCC, Clang, and MinGW Linux2Windows cross-compiler
+    - macOS 12 (Monterey) using Xcode
+    - Windows 11 using Visual Studio 2022 with MSVC and Clang
+
 #-------------------------------------------------------------------------------
 Author: Ling-Ti Kong, konglt(at)sjtu.edu.cn
-May 2020
+Aug 2021

tools/phonon/disp.cpp

@@ -1,10 +1,18 @@
-#include "string.h"
+
-#include "qnodes.h"
-#include "global.h"
 #include "phonon.h"
-#include "green.h"
+
-#include "timer.h"
+#include "dynmat.h"
+#include "global.h"
+#include "input.h"
 #include "kpath.h"
+#include "qnodes.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <string>
+#include <vector>
 
 /*------------------------------------------------------------------------------
  * Private method to evaluate the phonon dispersion curves
@@ -13,19 +21,22 @@ void Phonon::pdisp()
 {
   // ask the output file name and write the header.
    char str[MAXLINE];
-   for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+   puts("================================================================================");
+
 #ifdef UseSPG
   // ask method to generate q-lines
    int method = 2;
    printf("Please select your method to generate the phonon dispersion:\n");
    printf("  1. Manual, should always work;\n");
    printf("  2. Automatic, works only for 3D crystals (CMS49-299).\nYour choice [2]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f"));
+   input->read_stdin(str);
+   if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f"));
    method = 2 - method%2;
    printf("Your  selection: %d\n", method);
 #endif
    printf("\nPlease input the filename to output the dispersion data [pdisp.dat]:");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdisp.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str, "pdisp.dat");
    char *ptr = strtok(str," \t\n\r\f");
    char *fname = new char[strlen(ptr)+1];
    strcpy(fname,ptr);
@@ -45,9 +56,9 @@ void Phonon::pdisp()
       while (1){
          for (int i = 0; i < 3; ++i) qstr[i] = qend[i];
    
-         int quit = 0;
          printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
-         int n = count_words(fgets(str, MAXLINE, stdin));
+         input->read_stdin(str);
+         int n = count_words(str);
          ptr = strtok(str, " \t\n\r\f");
          if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
          else if (n >= 3){
@@ -56,14 +67,18 @@ void Phonon::pdisp()
            qstr[2] = atof(strtok(NULL, " \t\n\r\f"));
          }
      
-         do printf("Please input the end q-point in unit of B1->B3: ");
+         while ( 1 ){
-         while (count_words(fgets(str, MAXLINE, stdin)) < 3);
+            printf("Please input the end q-point in unit of B1->B3: ");
+            input->read_stdin(str);
+            if (count_words(str) >= 3) break;
+         }
          qend[0] = atof(strtok(str,  " \t\n\r\f"));
          qend[1] = atof(strtok(NULL, " \t\n\r\f"));
          qend[2] = atof(strtok(NULL, " \t\n\r\f"));
      
          printf("Please input the # of points along the line [%d]: ", nq);
-         if (count_words(fgets(str, MAXLINE, stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
+         input->read_stdin(str);
+         if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f"));
          nq = MAX(nq,2);
      
          double *qtmp = new double [3];
@@ -147,7 +162,7 @@ void Phonon::pdisp()
      printf("\nPhonon dispersion data are written to: %s, you can visualize the results\n", fname);
      printf("by invoking: `gnuplot pdisp.gnuplot; gv pdisp.eps`\n");
    }
-   for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+   puts("================================================================================");
  
    delete []fname;
    delete qnodes;

tools/phonon/dynmat.cpp

@@ -1,7 +1,16 @@
+
 #include "dynmat.h"
-#include "math.h"
+
-#include "version.h"
 #include "global.h"
+#include "input.h"
+#include "interpolate.h"
+#include "memory.h"
+#include "version.h"
+#include "zheevd.h"
+
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
 
 /* ----------------------------------------------------------------------------
  * Class DynMat stores the Dynamic Matrix read from the binary file from
@@ -9,6 +18,7 @@
  * ---------------------------------------------------------------------------- */
 DynMat::DynMat(int narg, char **arg)
 {
+   input = NULL;
    attyp = NULL;
    memory = NULL;
    M_inv_sqrt = NULL;
@@ -19,6 +29,8 @@ DynMat::DynMat(int narg, char **arg)
    attyp = NULL;
    basis = NULL;
    flag_reset_gamma = flag_skip = 0;
+   symprec = -1.;
+   int flag_save = 0;
 
    // analyze the command line options
    int iarg = 1;
@@ -29,9 +41,16 @@ DynMat::DynMat(int narg, char **arg)
       } else if (strcmp(arg[iarg], "-r") == 0){
          flag_reset_gamma = 1;
   
+      } else if (strcmp(arg[iarg], "-p") == 0){
+         if (++iarg >= narg) help();
+         else symprec = fabs(atof(arg[iarg]));
+  
       } else if (strcmp(arg[iarg], "-h") == 0){
          help();
   
+      } else if (strcmp(arg[iarg], "-save") == 0){
+        flag_save = 1;
+  
       } else {
          if (binfile) delete []binfile;
          int n = strlen(arg[iarg]) + 1;
@@ -43,6 +62,8 @@ DynMat::DynMat(int narg, char **arg)
    }
  
    ShowVersion();
+   input = new UserInput(flag_save);
+
    // get the binary file name from user input if not found in command line
    char str[MAXLINE];
    if (binfile == NULL) {
@@ -50,7 +71,7 @@ DynMat::DynMat(int narg, char **arg)
       printf("\n");
       do {
          printf("Please input the binary file name from fix_phonon: ");
-         fgets(str,MAXLINE,stdin);
+         input->read_stdin(str);
          ptr = strtok(str, " \n\t\r\f");
       } while (ptr == NULL);
   
@@ -137,17 +158,17 @@ DynMat::DynMat(int narg, char **arg)
       fclose(fp);
       exit(3);
    }
-   if (fread(basis[0], sizeof(double), fftdim, fp) != fftdim){
+   if (fread(basis[0], sizeof(double), fftdim, fp) != (size_t)fftdim){
       printf("\nError while reading basis info from file: %s\n", binfile);
       fclose(fp);
       exit(3);
    }
-   if (fread(&attyp[0], sizeof(int), nucell, fp) != nucell){
+   if (fread(&attyp[0], sizeof(int), nucell, fp) != (size_t)nucell){
       printf("\nError while reading atom types from file: %s\n", binfile);
       fclose(fp);
       exit(3);
    }
-   if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != nucell){
+   if (fread(&M_inv_sqrt[0], sizeof(double), nucell, fp) != (size_t)nucell){
       printf("\nError while reading atomic masses from file: %s\n", binfile);
       fclose(fp);
       exit(3);
@@ -159,6 +180,7 @@ DynMat::DynMat(int narg, char **arg)
 
    // initialize interpolation
    interpolate = new Interpolate(nx,ny,nz,fftdim2,DM_all);
+   interpolate->input = input;
    if (flag_reset_gamma) interpolate->reset_gamma();
 
    // Enforcing Austic Sum Rule
@@ -217,7 +239,7 @@ void DynMat::writeDMq(double *q)
       printf("\n");
       while ( 1 ){
          printf("Please input the filename to output the DM at selected q: ");
-         fgets(str,MAXLINE,stdin);
+         input->read_stdin(str);
          ptr = strtok(str, " \r\t\n\f");
          if (ptr) break;
       }
@@ -264,9 +286,9 @@ void DynMat::writeDMq(double *q, const double qr, FILE *fp)
 int DynMat::geteigen(double *egv, int flag)
 {
    char jobz, uplo;
-   integer n, lda, lwork, lrwork, *iwork, liwork, info;
+   int n, lda, lwork, lrwork, *iwork, liwork, info;
    doublecomplex *work;
-   doublereal *w = &egv[0], *rwork;
+   double *w = &egv[0], *rwork;
 
    n = fftdim;
    if (flag) jobz = 'V';
@@ -348,7 +370,7 @@ void DynMat::EnforceASR()
    char str[MAXLINE];
    int nasr = 20;
    if (nucell <= 1) nasr = 1;
-   printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
+   printf("\n================================================================================");
 
    // compute and display eigenvalues of Phi at gamma before ASR
    if (nucell > 100){
@@ -356,7 +378,7 @@ void DynMat::EnforceASR()
       fflush(stdout);
    }
 
-   double egvs[fftdim];
+   double *egvs = new double[fftdim];
    for (int i = 0; i < fftdim; ++i)
    for (int j = 0; j < fftdim; ++j) DM_q[i][j] = DM_all[0][i*fftdim+j];
    geteigen(egvs, 0);
@@ -370,11 +392,11 @@ void DynMat::EnforceASR()
  
    // ask for iterations to enforce ASR
    printf("Please input the # of iterations to enforce ASR [%d]: ", nasr);
-   fgets(str,MAXLINE,stdin);
+   input->read_stdin(str);
    char *ptr = strtok(str," \t\n\r\f");
    if (ptr) nasr = atoi(ptr);
    if (nasr < 1){
-      for (int i=0; i<80; i++) printf("="); printf("\n");
+      puts("================================================================================");
       return;
    }
  
@@ -439,9 +461,8 @@ void DynMat::EnforceASR()
       if (i%10 == 9) printf("\n");
       if (i == 99){ printf("...... (%d more skiped)", fftdim-100); break;}
    }
-   printf("\n");
+   delete[] egvs;
-   for (int i = 0; i < 80; ++i) printf("="); printf("\n\n");
+   puts("\n================================================================================\n");
- 
    return;
 }
 
@@ -468,7 +489,7 @@ void DynMat::real2rec()
  
    for (int i = 0; i < 9; ++i) ibasevec[i] *= vol;
  
-   printf("\n"); for (int i = 0; i < 80; ++i) printf("=");
+   printf("\n================================================================================");
    printf("\nBasis vectors of the unit cell in real space:");
    for (int i = 0; i < sysdim; ++i){
       printf("\n     A%d: ", i+1);
@@ -479,8 +500,7 @@ void DynMat::real2rec()
       printf("\n     B%d: ", i+1);
       for (int j = 0; j < sysdim; ++j) printf("%8.4f ", ibasevec[i*3+j]);
    }
-   printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
+   puts("\n================================================================================");
- 
    return;
 }
 
@@ -500,16 +520,17 @@ void DynMat::GaussJordan(int n, double *Mat)
    indxc = new int[n];
    indxr = new int[n];
    ipiv  = new int[n];
- 
+
+   irow = icol = -1;
    for (i = 0; i < n; ++i) ipiv[i] = 0;
    for (i = 0; i < n; ++i){
-      big = 0.;
+      big = 0.0;
       for (j = 0; j < n; ++j){
          if (ipiv[j] != 1){
             for (k = 0; k < n; ++k){
                if (ipiv[k] == 0){
                   idr = j * n + k;
-                  nmjk = abs(Mat[idr]);
+                  nmjk = fabs(Mat[idr]);
                   if (nmjk >= big){
                      big  = nmjk;
                      irow = j;
@@ -602,6 +623,9 @@ void DynMat::help()
    printf("              will also inform the code to skip all q-points that is in the vicinity\n");
    printf("              of the gamma point when evaluating phonon DOS and/or phonon dispersion.\n\n");
    printf("              By default, this is not set; and not expected for uncharged systems.\n\n");
+   printf("  -p prec     To define the precision for symmetry identification with spglib.\n");
+   printf("              By default, 1.e-3.\n\n");
+   printf("  -save       To record user input in `script.inp`, facilitating scripting.\n\n");
    printf("  -h          To print out this help info.\n\n");
    printf("  file        To define the filename that carries the binary dynamical matrice generated\n");
    printf("              by fix-phonon. If not provided, the code will ask for it.\n");
@@ -680,13 +704,13 @@ void DynMat::Define_Conversion_Factor()
  * ---------------------------------------------------------------------------- */
 void DynMat::ShowInfo()
 {
-   printf("\n"); for (int i = 0; i < 80; ++i) printf("="); printf("\n");
+   puts("\n================================================================================");
    printf("Dynamical matrix is read from file: %s\n", binfile);
    printf("The system size in three dimension: %d x %d x %d\n", nx, ny, nz);
    printf("Number of atoms per unit cell     : %d\n", nucell);
    printf("System dimension                  : %d\n", sysdim);
    printf("Boltzmann constant in used units  : %g\n", boltz);
-   for (int i = 0; i < 80; ++i) printf("="); printf("\n");
+   puts("================================================================================");
    return;
 }
 /* --------------------------------------------------------------------*/

tools/phonon/dynmat.h

@@ -1,11 +1,9 @@
 #ifndef DYNMAT_H
 #define DYNMAT_H
 
-#include "stdio.h"
+#include "zheevd.h"
-#include "stdlib.h"
+
-#include "string.h"
+#include <cstdio>
-#include "memory.h"
-#include "interpolate.h"
 
 class DynMat {
 public:
@@ -15,7 +13,7 @@ public:
 
   int nx, ny, nz, nucell;
   int sysdim, fftdim;
-  double eml2f, eml2fc;
+  double eml2f, eml2fc, symprec;
   char *funit;
 
   void getDMq(double *);
@@ -34,18 +32,18 @@ public:
   double **basis;
   int *attyp;
 
+  class UserInput *input;
+
 private:
 
   int flag_skip, flag_reset_gamma;
-  Interpolate *interpolate;
+  class Interpolate *interpolate;
-  
+  class Memory *memory;
-  Memory *memory;
 
-  int nasr;
   void EnforceASR();
 
   char *binfile, *dmfile;
-  double boltz, q[3];
+  double boltz;
 
   doublecomplex **DM_all;
 

tools/phonon/global.h

@@ -1,7 +1,7 @@
 #ifndef GLOBAL_H
 #define GLOBAL_H
 
-#define ZERO 1.e-8
+#define ZERO 1.0e-8
 #define MAXLINE 512
 
 #define MIN(a,b) ((a)>(b)?(b):(a))

tools/phonon/green.cpp

@@ -1,10 +1,10 @@
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
 #include "green.h"
+
+#include "memory.h"
+
 #include <complex>
-#include "global.h"
+#include <cmath>
+#include <cstdio>
 
 /*******************************************************************************
  * The class of Green is designed to evaluate the LDOS via the Green's Function
@@ -59,7 +59,6 @@ Green::Green(const int ntm, const int sdim, const int niter, const double min, c
    dw = (wmax - wmin)/double(nw-1);
    memory->create(alpha, sysdim,nit,  "Green_Green:alpha");
    memory->create(beta,  sysdim,nit+1,"Green_Green:beta");
-   //memory->create(ldos,  nw,sysdim, "Green_Green:ldos");
  
    // use Lanczos algorithm to diagonalize the Hessian
    Lanczos();
@@ -224,8 +223,6 @@ void Green::recursion()
 {
    // local variables
    std::complex<double> Z, rec_x, rec_x_inv;
-   std::complex<double> cunit = std::complex<double>(0.,1.);
- 
    double w = wmin;
  
    for (int i = 0; i < nw; ++i){

tools/phonon/green.h

@@ -1,8 +1,6 @@
 #ifndef GREEN_H
 #define GREEN_H
 
-#include "memory.h"
-
 class Green{
 public:
    Green(const int, const int, const int, const double, const double,
@@ -14,12 +12,11 @@ private:
    void Recursion();
    void recursion();
  
-   int ndos;
    double **ldos;
  
    int natom, iatom, sysdim, nit, nw, ndim;
    double dw, wmin, wmax, epson;
    double **alpha, **beta, **H; 
-   Memory *memory;
+   class Memory *memory;
 };
 #endif

tools/phonon/input.cpp

@@ -0,0 +1,35 @@
+#include "input.h"
+
+#include "global.h"
+
+/* -------------------------------------------------------------------
+ * Constructor. If flag = 1, output user inputs as script.inp
+ * ---------------------------------------------------------------- */ 
+UserInput::UserInput(int flag)
+{
+   fp = NULL;
+   if (flag) fp = fopen("script.inp", "w");
+
+   return;
+}
+
+/* -------------------------------------------------------------------
+ * Deconstructor. Output user inputs as required and clear workspace.
+ * ---------------------------------------------------------------- */ 
+UserInput::~UserInput()
+{
+   if (fp) fclose(fp);
+   fp = NULL;
+}
+
+/* -------------------------------------------------------------------
+ * Read stdin and keep a record of it.
+ * ---------------------------------------------------------------- */ 
+void UserInput::read_stdin(char *str)
+{
+   fgets(str, MAXLINE, stdin);
+   if (fp) fprintf(fp, "%s", str);
+
+   return;
+}
+/* ---------------------------------------------------------------- */

tools/phonon/input.h

@@ -0,0 +1,17 @@
+#ifndef INPUT_H
+#define INPUT_H
+
+#include <cstdio>
+
+class UserInput {
+public:
+   UserInput(int);
+   ~UserInput();
+
+   void read_stdin(char *);
+
+private:
+   FILE *fp;
+
+};
+#endif

tools/phonon/interpolate.cpp

@@ -1,6 +1,14 @@
+
 #include "interpolate.h"
-#include "math.h"
+
 #include "global.h"
+#include "input.h"
+#include "memory.h"
+#include "tricubic.h"
+
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
 /* ----------------------------------------------------------------------------
  * Constructor used to get info from caller, and prepare other necessary data
@@ -19,6 +27,7 @@ Interpolate::Interpolate(int nx, int ny, int nz, int ndm, doublecomplex **DM)
    data = DM;
    Dfdx = Dfdy = Dfdz = D2fdxdy = D2fdxdz = D2fdydz = D3fdxdydz = NULL;
    flag_reset_gamma = flag_allocated_dfs = 0;
+   input = NULL;
  
    return;
 }
@@ -265,17 +274,19 @@ void Interpolate::set_method()
 {
    char str[MAXLINE];
    int im = 1;
-   printf("\n");for(int i=0; i<80; i++) printf("=");
+   if (input == NULL) input = new UserInput(0);
-   printf("\nWhich interpolation method would you like to use?\n");
+
+   puts("\n================================================================================");
+   printf("Which interpolation method would you like to use?\n");
    printf("  1. Tricubic;\n  2. Trilinear;\n");
    printf("Your choice [1]: ");
-   fgets(str,MAXLINE,stdin);
+   input->read_stdin(str);
    char *ptr = strtok(str," \t\n\r\f");
    if (ptr) im = atoi(ptr);
  
    which =2-im%2;
    printf("Your  selection: %d\n", which);
-   for(int i=0; i<80; i++) printf("="); printf("\n\n");
+   puts("================================================================================\n");
  
    if (which == 1) tricubic_init();
  

tools/phonon/interpolate.h

@@ -1,16 +1,7 @@
 #ifndef INTERPOLATION_H
 #define INTERPOLATION_H
 
-#include "stdio.h"
+#include "zheevd.h"
-#include "stdlib.h"
-#include "string.h"
-#include "memory.h"
-#include "tricubic.h"
-
-extern "C"{
-#include "f2c.h"
-#include "clapack.h"
-}
 
 class Interpolate{
 public:
@@ -23,11 +14,13 @@ public:
  
    int UseGamma;
 
+   class UserInput *input;
+
 private:
    void tricubic_init();
    void tricubic(double *, doublecomplex *);
    void trilinear(double *, doublecomplex *);
-   Memory *memory;
+   class Memory *memory;
  
    int which;
    int Nx, Ny, Nz, Npt, ndim;

tools/phonon/kpath.cpp

@@ -1,11 +1,19 @@
-#include "global.h"
+
 #include "kpath.h"
 
+#include "global.h"
+#include "dynmat.h"
+#include "memory.h"
+#include "qnodes.h"
+
 #ifdef UseSPG
 extern "C"{
 #include "spglib.h"
 }
-#include "math.h"
+#include <cmath>
+#include <cstdio>
+#include <string>
+#include <vector>
 
 /* ----------------------------------------------------------------------------
  * Class kPath will help to find the high symmetry k-path for a given lattice.
@@ -47,11 +55,15 @@ kPath::kPath(DynMat *dm, QNodes *qn)
    for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim];
  
    // get the space group number
-   double symprec = 1.e-4, pos[num_atom][3];
+   double symprec = 1.0e-3;
+   double **pos;
+   memory->create(pos,num_atom,3,"kpath:pos");
+   if (dynmat->symprec > 0.0) symprec = dynmat->symprec;
+
    for (int i = 0; i < num_atom; ++i)
    for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j];
-   spgnum  = spg_get_international(symbol, latvec, pos, attyp, num_atom, symprec);
+   spgnum  = spg_get_international(symbol, latvec, (double (*)[3])pos, attyp, num_atom, symprec);
- 
+   memory->destroy(pos);
    return;
 }
 
@@ -61,7 +73,7 @@ kPath::kPath(DynMat *dm, QNodes *qn)
 void kPath::show_info()
 {
    // display the unit cell info read
-   for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n");
+   puts("--------------------------------------------------------------------------------");
    printf("The basis vectors of the unit cell:\n");
    for (int idim = 0; idim < 3; ++idim){
       printf("  A%d =", idim+1);
@@ -76,12 +88,10 @@ void kPath::show_info()
    if (num_atom > NUMATOM) printf("  ... (%d atoms omitted.)\n", num_atom-NUMATOM);
  
    printf("The space group number of your unit cell is: %d => %s\n", spgnum, symbol);
-   for (int ii = 0; ii < 80; ++ii) printf("-"); printf("\n");
+   puts("--------------------------------------------------------------------------------");
- 
    return;
 }
 
-
 /* ----------------------------------------------------------------------------
  * Free the memeory used by kPath.
  * ---------------------------------------------------------------------------- */
@@ -2765,9 +2775,16 @@ void kPath::show_path()
    if (q == NULL) return;
    int nbin = q->ndstr.size();
    if (nbin > 0){
-      printf("\nk-path for the current lattice will be:\n\t%s", q->ndstr[0].c_str());
+      puts("\n--------------------------------------------------------------------------------");
+      printf("k-path for the current lattice will be:\n  %s", q->ndstr[0].c_str());
       for (int is = 1; is < nbin; ++is) printf("-%s", q->ndstr[is].c_str());
-      printf("\n");
+
+      printf("\n\nThe fractional coordinates of these paths are:\n");
+      for (int is = 0; is < nbin-1; ++is)
+         printf("  [%6.4f %6.4f %6.4f] --> [%6.4f %6.4f %6.4f] (%s - %s)\n", q->qs[is][0],
+                q->qs[is][1], q->qs[is][2], q->qe[is][0], q->qe[is][1], q->qe[is][2],
+                q->ndstr[is].c_str(), q->ndstr[is+1].c_str() );
+      puts("--------------------------------------------------------------------------------");
    }
 
    return;

tools/phonon/kpath.h

@@ -4,14 +4,9 @@
 #ifndef KPATH_H
 #define KPATH_H
 
-#include "qnodes.h"
-#include "dynmat.h"
-#include "memory.h"
-
 class kPath{
 public:
-
+  kPath(class DynMat *, class QNodes *);
-  kPath(DynMat *, QNodes *);
   ~kPath();
 
   void kpath();
@@ -19,13 +14,11 @@ public:
   void show_info();
 
 private:
-
+  class Memory *memory;
-  Memory *memory;
+  class DynMat *dynmat;
-
+  class QNodes *q;
-  DynMat *dynmat;
-  QNodes *q;
   char symbol[11];
-  int spgnum, sysdim, fftdim, num_atom, *attyp;
+  int spgnum, sysdim, num_atom, *attyp;
   double latvec[3][3], **atpos;
 
 };

tools/phonon/main.cpp

@@ -1,10 +1,6 @@
-#include "stdio.h"
-#include "stdlib.h"
 #include "dynmat.h"
 #include "phonon.h"
 
-using namespace std;
-
 int main(int argc, char** argv)
 {
 

tools/phonon/memory.cpp

@@ -1,8 +1,8 @@
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
 #include "memory.h"
 
+#include <cstdio>
+#include <cstdlib>
+
 /* ----------------------------------------------------------------------
    safe malloc 
 ------------------------------------------------------------------------- */

tools/phonon/memory.h

@@ -4,18 +4,16 @@
 #define __STDC_LIMIT_MACROS
 #define __STDC_FORMAT_MACROS
 
-#include "stdio.h"
+#include <cstdlib>
-#include "stdlib.h"
+#include <stdint.h>
-#include "limits.h"
+#include <inttypes.h>
-#include "stdint.h"
-#include "inttypes.h"
 
 typedef int64_t bigint;
 #define BIGINT_FORMAT "%" PRId64
 #define ATOBIGINT atoll
 
 class Memory {
- public:
+public:
   Memory(){};
 
   void *smalloc(bigint n, const char *);
@@ -24,11 +22,11 @@ class Memory {
   void fail(const char *);
 
 /* ----------------------------------------------------------------------
-   create a 1d array 
+   create a 1d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *create(TYPE *&array, int n, const char *name)
+  TYPE *create(TYPE *&array, int n, const char *name)
     {
       bigint nbytes = sizeof(TYPE) * n;
       array = (TYPE *) smalloc(nbytes,name);
@@ -36,42 +34,42 @@ class Memory {
     };
 
   template <typename TYPE>
-    TYPE **create(TYPE **&array, int n, const char *name) {fail(name);}
+  TYPE **create(TYPE **&, int, const char *name) {fail(name);}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *grow(TYPE *&array, int n, const char *name) 
+  TYPE *grow(TYPE *&array, int n, const char *name)
     {
       if (array == NULL) return create(array,n,name);
-      
+
       bigint nbytes = sizeof(TYPE) * n;
       array = (TYPE *) srealloc(array,nbytes,name);
       return array;
     };
 
   template <typename TYPE>
-    TYPE **grow(TYPE **&array, int n, const char *name) {fail(name);}
+  TYPE **grow(TYPE **&, int, const char *name) {fail(name);}
 
 /* ----------------------------------------------------------------------
-   destroy a 1d array 
+   destroy a 1d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE *array) 
+  void destroy(TYPE *array)
     {
       sfree(array);
     };
 
 /* ----------------------------------------------------------------------
-   create a 1d array with index from nlo to nhi inclusive 
+   create a 1d array with index from nlo to nhi inclusive
    cannot grow it
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name) 
+  TYPE *create1d_offset(TYPE *&array, int nlo, int nhi, const char *name)
     {
       bigint nbytes = sizeof(TYPE) * (nhi-nlo+1);
       array = (TYPE *) smalloc(nbytes,name);
@@ -80,76 +78,76 @@ class Memory {
     }
 
   template <typename TYPE>
-    TYPE **create1d_offset(TYPE **&array, int nlo, int nhi, const char *name)
+  TYPE **create1d_offset(TYPE **&, int, int, const char *name)
     {fail(name);}
 
 /* ----------------------------------------------------------------------
-   destroy a 1d array with index offset 
+   destroy a 1d array with index offset
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy1d_offset(TYPE *array, int offset) 
+  void destroy1d_offset(TYPE *array, int offset)
     {
       if (array) sfree(&array[offset]);
     }
 
 /* ----------------------------------------------------------------------
-   create a 2d array 
+   create a 2d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **create(TYPE **&array, int n1, int n2, const char *name) 
+  TYPE **create(TYPE **&array, int n1, int n2, const char *name)
     {
       bigint nbytes = sizeof(TYPE) * n1*n2;
       TYPE *data = (TYPE *) smalloc(nbytes,name);
       nbytes = sizeof(TYPE *) * n1;
       array = (TYPE **) smalloc(nbytes,name);
-      
+
       int n = 0;
       for (int i = 0; i < n1; i++) {
-	array[i] = &data[n];
+        array[i] = &data[n];
-	n += n2;
+        n += n2;
       }
       return array;
     }
 
   template <typename TYPE>
-    TYPE ***create(TYPE ***&array, int n1, int n2, const char *name)
+  TYPE ***create(TYPE ***&, int, int, const char *name)
     {fail(name);}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1st dim of a 2d array
    last dim must stay the same
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **grow(TYPE **&array, int n1, int n2, const char *name) 
+  TYPE **grow(TYPE **&array, int n1, int n2, const char *name)
     {
       if (array == NULL) return create(array,n1,n2,name);
-      
+
       bigint nbytes = sizeof(TYPE) * n1*n2;
       TYPE *data = (TYPE *) srealloc(array[0],nbytes,name);
       nbytes = sizeof(TYPE *) * n1;
       array = (TYPE **) srealloc(array,nbytes,name);
-      
+
       int n = 0;
       for (int i = 0; i < n1; i++) {
-	array[i] = &data[n];
+        array[i] = &data[n];
-	n += n2;
+        n += n2;
       }
       return array;
     }
 
   template <typename TYPE>
-    TYPE ***grow(TYPE ***&array, int n1, int n2, const char *name)
+  TYPE ***grow(TYPE ***&, int, int, const char *name)
     {fail(name);}
 
 /* ----------------------------------------------------------------------
-   destroy a 2d array 
+   destroy a 2d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE **array)
+  void destroy(TYPE **array)
     {
       if (array == NULL) return;
       sfree(array[0]);
@@ -157,42 +155,11 @@ class Memory {
     }
 
 /* ----------------------------------------------------------------------
-   create a 2d array with 2nd index from n2lo to n2hi inclusive 
+   create a 3d array
-   cannot grow it
+   ------------------------------------------------------------------------- */
-------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi,
+  TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name)
-			   const char *name)
-    {
-      int n2 = n2hi - n2lo + 1;
-      create(array,n1,n2,name);
-      for (int i = 0; i < n1; i++) array[i] -= n2lo;
-      return array;
-    }
-
-  template <typename TYPE>
-    TYPE ***create2d_offset(TYPE ***&array, int n1, int n2lo, int n2hi,
-			    const char *name) {fail(name);}
-
-/* ----------------------------------------------------------------------
-   destroy a 2d array with 2nd index offset 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    void destroy2d_offset(TYPE **array, int offset)
-    {
-      if (array == NULL) return;
-      sfree(&array[0][offset]);
-      sfree(array);
-    }
-
-/* ----------------------------------------------------------------------
-   create a 3d array 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) 
     {
       bigint nbytes = sizeof(TYPE) * n1*n2*n3;
       TYPE *data = (TYPE *) smalloc(nbytes,name);
@@ -200,62 +167,62 @@ class Memory {
       TYPE **plane = (TYPE **) smalloc(nbytes,name);
       nbytes = sizeof(TYPE **) * n1;
       array = (TYPE ***) smalloc(nbytes,name);
-      
+
       int i,j;
       int n = 0;
       for (i = 0; i < n1; i++) {
-	array[i] = &plane[i*n2];
+        array[i] = &plane[i*n2];
-	for (j = 0; j < n2; j++) {
+        for (j = 0; j < n2; j++) {
-	  plane[i*n2+j] = &data[n];
+          plane[i*n2+j] = &data[n];
-	  n += n3;
+          n += n3;
-	}
+        }
       }
       return array;
     }
 
   template <typename TYPE>
-    TYPE ****create(TYPE ****&array, int n1, int n2, int n3, const char *name) 
+  TYPE ****create(TYPE ****&, int, int, int, const char *name)
     {fail(name);}
 
 /* ----------------------------------------------------------------------
    grow or shrink 1st dim of a 3d array
    last 2 dims must stay the same
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name) 
+  TYPE ***grow(TYPE ***&array, int n1, int n2, int n3, const char *name)
     {
       if (array == NULL) return create(array,n1,n2,n3,name);
-      
+
       bigint nbytes = sizeof(TYPE) * n1*n2*n3;
       TYPE *data = (TYPE *) srealloc(array[0][0],nbytes,name);
       nbytes = sizeof(TYPE *) * n1*n2;
       TYPE **plane = (TYPE **) srealloc(array[0],nbytes,name);
       nbytes = sizeof(TYPE **) * n1;
       array = (TYPE ***) srealloc(array,nbytes,name);
-      
+
       int i,j;
       int n = 0;
       for (i = 0; i < n1; i++) {
-	array[i] = &plane[i*n2];
+        array[i] = &plane[i*n2];
-	for (j = 0; j < n2; j++) {
+        for (j = 0; j < n2; j++) {
-	  plane[i*n2+j] = &data[n];
+          plane[i*n2+j] = &data[n];
-	  n += n3;
+          n += n3;
-	}
+        }
       }
       return array;
     }
 
   template <typename TYPE>
-    TYPE ****grow(TYPE ****&array, int n1, int n2, int n3, const char *name) 
+  TYPE ****grow(TYPE ****&, int, int, int, const char *name)
     {fail(name);}
 
 /* ----------------------------------------------------------------------
-   destroy a 3d array 
+   destroy a 3d array
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    void destroy(TYPE ***array) 
+  void destroy(TYPE ***array)
     {
       if (array == NULL) return;
       sfree(array[0][0]);
@@ -263,168 +230,23 @@ class Memory {
       sfree(array);
     }
 
-/* ----------------------------------------------------------------------
-   create a 3d array with 1st index from n1lo to n1hi inclusive 
-   cannot grow it
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, 
-			    int n2, int n3, const char *name)
-    {
-      int n1 = n1hi - n1lo + 1;
-      create(array,n1,n2,n3,name);
-      array -= n1lo;
-      return array;
-    }
-
-  template <typename TYPE>
-    TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, 
-			     int n2, int n3, const char *name)
-    {fail(name);}
-
-/* ----------------------------------------------------------------------
-   free a 3d array with 1st index offset 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    void destroy3d_offset(TYPE ***array, int offset)
-    {
-      if (array) destroy(&array[offset]);
-    }
-
-/* ----------------------------------------------------------------------
-   create a 3d array with
-   1st index from n1lo to n1hi inclusive,
-   2nd index from n2lo to n2hi inclusive,
-   3rd index from n3lo to n3hi inclusive 
-   cannot grow it
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    TYPE ***create3d_offset(TYPE ***&array, int n1lo, int n1hi, 
-			    int n2lo, int n2hi, int n3lo, int n3hi,
-			    const char *name)
-    {
-      int n1 = n1hi - n1lo + 1;
-      int n2 = n2hi - n2lo + 1;
-      int n3 = n3hi - n3lo + 1;
-      create(array,n1,n2,n3,name);
-      
-      for (int i = 0; i < n1*n2; i++) array[0][i] -= n3lo;
-      for (int i = 0; i < n1; i++) array[i] -= n2lo;
-      array -= n1lo;
-      return array;
-    }
-
-  template <typename TYPE>
-    TYPE ****create3d_offset(TYPE ****&array, int n1lo, int n1hi, 
-			     int n2lo, int n2hi, int n3lo, int n3hi,
-			     const char *name)
-    {fail(name);}
-
-/* ----------------------------------------------------------------------
-   free a 3d array with all 3 indices offset 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    void destroy3d_offset(TYPE ***array, 
-			  int n1_offset, int n2_offset, int n3_offset)
-    {
-      if (array == NULL) return;
-      sfree(&array[n1_offset][n2_offset][n3_offset]);
-      sfree(&array[n1_offset][n2_offset]);
-      sfree(&array[n1_offset]);
-    }
-
-/* ----------------------------------------------------------------------
-   create a 4d array 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    TYPE ****create(TYPE ****&array, int n1, int n2, int n3, int n4,
-		    const char *name)
-    {
-      bigint nbytes = sizeof(TYPE) * n1*n2*n3*n4;
-      TYPE *data = (double *) smalloc(nbytes,name);
-      nbytes = sizeof(TYPE *) * n1*n2*n3;
-      TYPE **cube = (double **) smalloc(nbytes,name);
-      nbytes = sizeof(TYPE **) * n1*n2;
-      TYPE ***plane = (double ***) smalloc(nbytes,name);
-      nbytes = sizeof(TYPE ***) * n1;
-      array = (double ****) smalloc(nbytes,name);
-      
-      int i,j,k;
-      int n = 0;
-      for (i = 0; i < n1; i++) {
-	array[i] = &plane[i*n2];
-	for (j = 0; j < n2; j++) {
-	  plane[i*n2+j] = &cube[i*n2*n3+j*n3];
-	  for (k = 0; k < n3; k++) {
-	    cube[i*n2*n3+j*n3+k] = &data[n];
-	    n += n4;
-	  }
-	}
-      }
-      return array;
-    }
-
-  template <typename TYPE>
-    TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4,
-		     const char *name)
-    {fail(name);}
-
-/* ----------------------------------------------------------------------
-   destroy a 4d array 
-------------------------------------------------------------------------- */
-
-  template <typename TYPE>
-    void destroy(TYPE ****array)
-    {
-      if (array == NULL) return;
-      sfree(array[0][0][0]);
-      sfree(array[0][0]);
-      sfree(array[0]);
-      sfree(array);
-    }
-
 /* ----------------------------------------------------------------------
    memory usage of arrays, including pointers
-------------------------------------------------------------------------- */
+   ------------------------------------------------------------------------- */
 
   template <typename TYPE>
-    bigint usage(TYPE *array, int n)
+  bigint usage(TYPE *, int n)
     {
       bigint bytes = sizeof(TYPE) * n;
       return bytes;
     }
 
   template <typename TYPE>
-    bigint usage(TYPE **array, int n1, int n2)
+  bigint usage(TYPE **, int n1, int n2)
     {
       bigint bytes = sizeof(TYPE) * n1*n2;
       bytes += sizeof(TYPE *) * n1;
       return bytes;
     }
-
-  template <typename TYPE>
-    bigint usage(TYPE ***array, int n1, int n2, int n3)
-    {
-      bigint bytes = sizeof(TYPE) * n1*n2*n3;
-      bytes += sizeof(TYPE *) * n1*n2;
-      bytes += sizeof(TYPE **) * n1;
-      return bytes;
-    }
-
-  template <typename TYPE>
-    bigint usage(TYPE ****array, int n1, int n2, int n3, int n4)
-    {
-      bigint bytes = sizeof(TYPE) * n1*n2*n3*n4;
-      bytes += sizeof(TYPE *) * n1*n2*n3;
-      bytes += sizeof(TYPE **) * n1*n2;
-      bytes += sizeof(TYPE ***) * n1;
-      return bytes;
-    }
 };
-
 #endif

tools/phonon/phonon.cpp

@@ -1,9 +1,18 @@
-#include <vector>
+
-#include "string.h"
 #include "phonon.h"
-#include "green.h"
+
-#include "timer.h"
 #include "global.h"
+#include "dynmat.h"
+#include "green.h"
+#include "input.h"
+#include "memory.h"
+#include "timer.h"
+#include "zheevd.h"
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
 
 #ifdef UseSPG
 extern "C"{
@@ -27,6 +36,7 @@ Phonon::Phonon(DynMat *dm)
    dynmat = dm;
    sysdim = dynmat->sysdim;
    ndim   = dynmat->fftdim;
+   input  = dm->input;
    dos  = NULL;
    ldos = NULL;
    qpts = NULL;
@@ -42,10 +52,7 @@ Phonon::Phonon(DynMat *dm)
    // display the menu
    char str[MAXLINE];
    while ( 1 ){
-      printf("\n");
+      puts("\n===================================== Menu =====================================");
-      for (int i = 0; i < 37; ++i) printf("=");
-      printf(" Menu ");
-      for (int i = 0; i < 37; ++i) printf("="); printf("\n");
       printf("  1. Phonon DOS evaluation;\n");
       printf("  2. Phonon dispersion curves;\n");
       printf("  3. Dynamical matrix at arbitrary q;\n");
@@ -64,9 +71,10 @@ Phonon::Phonon(DynMat *dm)
       // read user choice
       int job = 0;
       printf("Your choice [0]: ");
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) job = atoi(strtok(str," \t\n\r\f"));
+      input->read_stdin(str);
+      if (count_words(str) > 0) job = atoi(strtok(str," \t\n\r\f"));
       printf("\nYour  selection: %d\n", job);
-      for (int i = 0; i < 80; ++i) printf("=");printf("\n\n");
+      puts("================================================================================\n");
   
       // now to do the job according to user's choice
       if      (job == 1) pdos();
@@ -138,7 +146,8 @@ void Phonon::pdos()
    // Now to ask for the output frequency range
    printf("\nThe frequency range of all q-points are: [%g %g]\n", fmin, fmax);
    printf("Please input the desired range to get DOS [%g %g]: ", fmin, fmax);
-   if (count_words(fgets(str,MAXLINE,stdin)) >= 2){
+   input->read_stdin(str);
+   if (count_words(str) >= 2){
       fmin = atof(strtok(str," \t\n\r\f"));
       fmax = atof(strtok(NULL," \t\n\r\f"));
    }
@@ -147,7 +156,8 @@ void Phonon::pdos()
  
    ndos = 201;
    printf("Please input the number of intervals [%d]: ", ndos);
-   if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
+   input->read_stdin(str);
+   if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
    ndos += (ndos+1)%2; ndos = MAX(2,ndos);
  
    
@@ -170,7 +180,8 @@ void Phonon::pdos()
  
    // smooth dos ?
    printf("Would you like to smooth the phonon dos? (y/n)[n]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) > 0){
+   input->read_stdin(str);
+   if (count_words(str) > 0){
       char *flag = strtok(str," \t\n\r\f");
       if (strcmp(flag,"y") == 0 || strcmp(flag,"Y") == 0) smooth(dos, ndos);
    }
@@ -194,7 +205,8 @@ void Phonon::writeDOS()
    char str[MAXLINE];
    // now to output the phonon DOS
    printf("\nPlease input the filename to write DOS [pdos.dat]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "pdos.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str, "pdos.dat");
    char *fname = strtok(str," \t\n\r\f");
  
    printf("The total phonon DOS will be written to file: %s\n", fname);
@@ -234,7 +246,7 @@ void Phonon::writeLDOS()
    const double one3 = 1./double(sysdim);
    char str[MAXLINE];
    for (int ilocal = 0; ilocal < nlocal; ++ilocal){
-      sprintf(str,"pldos_%d.dat", locals[ilocal]);
+      snprintf(str, MAXLINE-1, "pldos_%d.dat", locals[ilocal]);
       char *fname = strtok(str," \t\n\r\f");
   
       FILE *fp = fopen(fname, "w"); fname = NULL;
@@ -281,7 +293,7 @@ void Phonon::ldos_rsgf()
   
       fmin = fmax = egvs[0];
       for (int i = 1; i < ndim; ++i){fmin = MIN(fmin, egvs[i]); fmax = MAX(fmax, egvs[i]);}
-      delete []egvs;
+      delete[] egvs;
 
    } else {
      fmin = 0.; fmax = 20.;
@@ -297,7 +309,8 @@ void Phonon::ldos_rsgf()
       printf("\nThere are %d atoms in each unit cell of your lattice.\n", dynmat->nucell);
       printf("Please input the index/index range/index range and increment of atom(s)\n");
       printf("in the unit cell to evaluate LDOS, q to exit [%d]: ", ik);
-      int nr = count_words( fgets(str,MAXLINE,stdin) );
+      input->read_stdin(str);
+      int nr = count_words(str);
       if (nr < 1){
          istr = iend = ik;
          iinc = 1;
@@ -327,7 +340,8 @@ void Phonon::ldos_rsgf()
       }
   
       printf("Please input the frequency range to evaluate LDOS [%g %g]: ", fmin, fmax);
-      if (count_words(fgets(str,MAXLINE,stdin)) >= 2){
+      input->read_stdin(str);
+      if (count_words(str) >= 2){
          fmin = atof(strtok(str," \t\n\r\f"));
          fmax = atof(strtok(NULL," \t\n\r\f"));
       }
@@ -335,16 +349,19 @@ void Phonon::ldos_rsgf()
       printf("The frequency range for your LDOS is [%g %g].\n", fmin, fmax);
   
       printf("Please input the desired number of points in LDOS [%d]: ", ndos);
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
+      input->read_stdin(str);
+      if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
       if (ndos < 2) break;
       ndos += (ndos+1)%2;
   
       printf("Please input the maximum # of Lanczos iterations  [%d]: ", nit);
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) nit = atoi(strtok(str," \t\n\r\f"));
+      input->read_stdin(str);
+      if (count_words(str) > 0) nit = atoi(strtok(str," \t\n\r\f"));
       if (nit < 1) break;
   
       printf("Please input the value of epsilon for delta-function [%g]: ", eps);
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) eps = atof(strtok(str," \t\n\r\f"));
+      input->read_stdin(str);
+      if (count_words(str) > 0) eps = atof(strtok(str," \t\n\r\f"));
       if (eps <= 0.) break;
       
       // prepare array for local pdos
@@ -395,8 +412,11 @@ void Phonon::dmanyq()
 {
    char str[MAXLINE];
    double q[3];
-   do printf("Please input the q-point to output the dynamical matrix:");
+   while ( 1 ){
-   while (count_words(fgets(str,MAXLINE,stdin)) < 3);
+      printf("Please input the q-point to output the dynamical matrix: ");
+      input->read_stdin(str);
+      if (count_words(str) >= 3) break;
+   }
    q[0] = atof(strtok(str," \t\n\r\f"));
    q[1] = atof(strtok(NULL," \t\n\r\f"));
    q[2] = atof(strtok(NULL," \t\n\r\f"));
@@ -413,11 +433,13 @@ void Phonon::dmanyq()
 void Phonon::vfanyq()
 {
    char str[MAXLINE];
-   double q[3], egvs[ndim];
+   double q[3];
+   double *egvs = new double[ndim];
    
    while ( 1 ){
       printf("Please input the q-point to compute the frequencies, q to exit: ");
-      if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
+      input->read_stdin(str);
+      if (count_words(str) < 3) break;
   
       q[0] = atof(strtok(str, " \t\n\r\f"));
       q[1] = atof(strtok(NULL," \t\n\r\f"));
@@ -427,9 +449,11 @@ void Phonon::vfanyq()
       dynmat->geteigen(egvs, 0);
       printf("q-point: [%lg %lg %lg], ", q[0], q[1], q[2]);
       printf("vibrational frequencies at this q-point:\n");
-      for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]); printf("\n\n");
+      for (int i = 0; i < ndim; ++i) printf("%lg ", egvs[i]);
+      printf("\n\n");
    }
- 
+
+   delete[] egvs;
    return;
 }
 
@@ -439,15 +463,18 @@ void Phonon::vfanyq()
 void Phonon::vecanyq()
 {
    char str[MAXLINE];
-   double q[3], egvs[ndim];
+   double q[3];
+   double *egvs = new double[ndim];
    doublecomplex **eigvec = dynmat->DM_q;
    printf("Please input the filename to output the result [eigvec.dat]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str,"eigvec.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str,"eigvec.dat");
    FILE *fp = fopen(strtok(str," \t\n\r\f"), "w");
  
    while ( 1 ){
       printf("Please input the q-point to compute the frequencies, q to exit: ");
-      if (count_words(fgets(str,MAXLINE,stdin)) < 3) break;
+      input->read_stdin(str);
+      if (count_words(str) < 3) break;
   
       q[0] = atof(strtok(str, " \t\n\r\f"));
       q[1] = atof(strtok(NULL," \t\n\r\f"));
@@ -475,6 +502,7 @@ void Phonon::vecanyq()
       fprintf(fp,"\n");
    }
    fclose(fp);
+   delete[] egvs;
    eigvec = NULL;
    return;
 }
@@ -488,7 +516,8 @@ void Phonon::DMdisp()
    char str[MAXLINE];
  
    printf("Please input the filename to output the DM data [DMDisp.dat]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "DMDisp.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str, "DMDisp.dat");
    char *fname = strtok(str," \t\n\r\f");
  
    FILE *fp = fopen(fname, "w"); fname = NULL;
@@ -503,7 +532,8 @@ void Phonon::DMdisp()
       for (int i = 0; i < 3; ++i) qstr[i] = qend[i];
   
       printf("\nPlease input the start q-point in unit of B1->B3, q to exit [%g %g %g]: ", qstr[0], qstr[1], qstr[2]);
-      int n = count_words(fgets(str,MAXLINE,stdin));
+      input->read_stdin(str);
+      int n = count_words(str);
       char *ptr = strtok(str," \t\n\r\f");
       if ((n == 1) && (strcmp(ptr,"q") == 0)) break;
       else if (n >= 3){
@@ -512,14 +542,18 @@ void Phonon::DMdisp()
          qstr[2] = atof(strtok(NULL," \t\n\r\f"));
       }
   
-      do printf("Please input the end q-point in unit of B1->B3: ");
+      while ( 1 ){
-      while (count_words(fgets(str,MAXLINE,stdin)) < 3);
+         printf("Please input the end q-point in unit of B1->B3: ");
+         input->read_stdin(str);
+         if (count_words(str) >= 3) break;
+      }
       qend[0] = atof(strtok(str," \t\n\r\f"));
       qend[1] = atof(strtok(NULL," \t\n\r\f"));
       qend[2] = atof(strtok(NULL," \t\n\r\f"));
   
       printf("Please input the # of points along the line [%d]: ", nq);
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0) nq = atoi(strtok(str," \t\n\r\f"));
+      input->read_stdin(str);
+      if (count_words(str) > 0) nq = atoi(strtok(str," \t\n\r\f"));
       nq = MAX(nq,2);
   
       for (int i=0; i<3; i++) qinc[i] = (qend[i]-qstr[i])/double(nq-1);
@@ -588,7 +622,8 @@ void Phonon::therm()
    char str[MAXLINE];
  
    printf("\nPlease input the filename to output thermal properties [therm.dat]:");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "therm.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str, "therm.dat");
    char *fname = strtok(str," \t\n\r\f");
    FILE *fp = fopen(fname, "a"); fname = NULL;
    // header line 
@@ -630,7 +665,8 @@ void Phonon::therm()
       fprintf(fp,"%lg %lg %lg %lg %lg %lg\n", T, Uvib, Svib, Fvib, ZPE, Cvib);
   
       printf("Please input the desired temperature (K), enter to exit: ");
-      if (count_words(fgets(str,MAXLINE,stdin)) < 1) break;
+      input->read_stdin(str);
+      if (count_words(str) < 1) break;
       T = atof(strtok(str," \t\n\r\f"));
   
    } while (T > 0.);
@@ -646,12 +682,14 @@ void Phonon::local_therm()
 {
    char str[MAXLINE];
    printf("\nWould you like to compute the local thermal properties (y/n)[n]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) return;
+   input->read_stdin(str);
+   if (count_words(str) < 1) return;
    char *ptr = strtok(str," \t\n\r\f");
    if (strcmp(ptr,"y") != 0 && strcmp(ptr, "Y") != 0 && strcmp(ptr, "yes") != 0) return;
  
    printf("Please input the filename to output vibrational thermal info [localtherm.dat]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) < 1) strcpy(str, "localtherm.dat");
+   input->read_stdin(str);
+   if (count_words(str) < 1) strcpy(str, "localtherm.dat");
  
    FILE *fp = fopen(strtok(str," \t\n\r\f"), "w");
    fprintf(fp,"# atom Temp  U_vib (eV)    S_vib (kB)    F_vib (eV)    C_vib (kB)     ZPE (eV)\n");
@@ -672,7 +710,8 @@ void Phonon::local_therm()
    while ( 1 ){
       printf("\nPlease input the temperature at which to evaluate the local vibrational\n");
       printf("thermal properties, non-positive number to exit [%g]: ", T);
-      if (count_words(fgets(str,MAXLINE,stdin)) > 0){
+      input->read_stdin(str);
+      if (count_words(str) > 0){
          T = atoi(strtok(str," \t\n\r\f"));
          if (T <= 0.) break;
       }
@@ -765,7 +804,8 @@ void Phonon::QMesh()
    printf("\nThe q-mesh size from the read dynamical matrix is: %d x %d x %d\n", nx, ny, nz);
    printf("A denser mesh can be interpolated, but NOTE a too dense mesh can cause segmentation fault.\n");
    printf("Please input your desired q-mesh size [%d %d %d]: ", nx, ny, nz);
-   if (count_words(fgets(str,MAXLINE,stdin)) >= 3){
+   input->read_stdin(str);
+   if (count_words(str) >= 3){
       nx = atoi(strtok(str," \t\n\r\f"));
       ny = atoi(strtok(NULL," \t\n\r\f"));
       nz = atoi(strtok(NULL," \t\n\r\f"));
@@ -780,7 +820,8 @@ void Phonon::QMesh()
    int method = 2;
    printf("Please select your method to generate the q-points:\n");
    printf("  1. uniform;\n  2. Monkhost-Pack mesh;\nYour choice [2]: ");
-   if (count_words(fgets(str,MAXLINE,stdin)) > 0) method = atoi(strtok(str," \t\n\r\f"));
+   input->read_stdin(str);
+   if (count_words(str) > 0) method = atoi(strtok(str," \t\n\r\f"));
    method = 2 - method%2;
    printf("Your  selection: %d\n", method);
 #endif
@@ -831,7 +872,7 @@ void Phonon::QMesh()
       for (int idim = 0; idim < sysdim; ++idim) atpos[i][idim] = dynmat->basis[i][idim];
   
       // display the unit cell info read
-      printf("\n");for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+      puts("\n================================================================================");
       printf("The basis vectors of the unit cell:\n");
       for (int idim = 0; idim < 3; ++idim)
          printf("  A%d = %lg %lg %lg\n", idim+1, latvec[0][idim], latvec[1][idim], latvec[2][idim]);
@@ -845,14 +886,21 @@ void Phonon::QMesh()
       mesh[0] = nx; mesh[1] = ny; mesh[2] = nz;
       shift[0] = shift[1] = shift[2] = 0;
       int num_grid = mesh[0]*mesh[1]*mesh[2];
-      int grid_point[num_grid][3], map[num_grid];
+      int **grid_point;
-      double symprec = 1.e-4, pos[num_atom][3];
+      memory->create(grid_point, num_grid, 3, "phonon:grid_point");
+      int *map = new int[num_grid];
+      double symprec = 1.0e-3;
+      double **pos;
+      memory->create(pos, num_atom, 3, "phonon:pos");
+      if (dynmat->symprec > 0.) symprec = dynmat->symprec;
   
       for (int i = 0; i < num_atom; ++i)
-      for (int j = 0; j < 3; ++j) pos[i][j] = atpos[i][j];
+          for (int j = 0; j < 3; ++j)
+              pos[i][j] = atpos[i][j];
   
       // if spglib >= 1.0.3 is used
-      nq = spg_get_ir_reciprocal_mesh(grid_point, map, mesh, shift, is_time_reversal, latvec, pos, attyp, num_atom, symprec);
+      nq = spg_get_ir_reciprocal_mesh((int (*)[3])grid_point, map, mesh, shift, is_time_reversal,
+                                      latvec, (double (*)[3])pos, attyp, num_atom, symprec);
   
       memory->create(wt,   nq, "QMesh:wt");
       memory->create(qpts, nq,3,"QMesh:qpts");
@@ -873,11 +921,14 @@ void Phonon::QMesh()
             qpts[numq][2] = double(grid_point[i][2])/double(mesh[2]);
             numq++;
          }
-         wt[iq2idx[iq]] += 1.;
+         wt[iq2idx[iq]] += 1.0;
       }
-      delete []iq2idx;
+      delete[] iq2idx;
-  
+      delete[] map;
-      double wsum = 0.;
+      memory->destroy(grid_point);
+      memory->destroy(pos);
+
+      double wsum = 0.0;
       for (int iq = 0; iq < nq; ++iq) wsum += wt[iq];
       for (int iq = 0; iq < nq; ++iq) wt[iq] /= wsum;
       
@@ -898,7 +949,8 @@ void Phonon::ldos_egv()
    char str[MAXLINE], *ptr;
    printf("\nThe # of atoms per cell is: %d, please input the atom IDs to compute\n", dynmat->nucell);
    printf("local PDOS, IDs begin with 0: ");
-   int nmax = count_words(fgets(str,MAXLINE,stdin));
+   input->read_stdin(str);
+   int nmax = count_words(str);
    if (nmax < 1) return;
  
    memory->destroy(locals);
@@ -920,7 +972,8 @@ void Phonon::ldos_egv()
  
    fmin = 0.; fmax = 10.;
    printf("Please input the freqency (nv, THz) range to compute PDOS [%g %g]: ", fmin, fmax);
-   if (count_words(fgets(str,MAXLINE,stdin)) >= 2) {
+   input->read_stdin(str);
+   if (count_words(str) >= 2) {
       fmin = atof(strtok(str," \t\n\r\f"));
       fmax = atof(strtok(NULL," \t\n\r\f"));
    }
@@ -928,7 +981,8 @@ void Phonon::ldos_egv()
  
    ndos = 201;
    printf("Please input your desired # of points in PDOS [%d]: ", ndos);
-   if (count_words(fgets(str,MAXLINE,stdin)) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
+   input->read_stdin(str);
+   if (count_words(str) > 0) ndos = atoi(strtok(str," \t\n\r\f"));
    if (ndos < 2) return;
    ndos += (ndos+1)%2;
  
@@ -957,7 +1011,8 @@ void Phonon::ldos_egv()
    Timer *time = new Timer();
  
    // memory and pointer for eigenvalues and eigenvectors
-   double egval[ndim], offset=fmin-0.5*df;
+   double offset=fmin-0.5*df;
+   double *egval = new double[ndim];
    doublecomplex **egvec = dynmat->DM_q;
  
    printf("\nNow to compute the phonons and DOSs "); fflush(stdout);
@@ -985,6 +1040,7 @@ void Phonon::ldos_egv()
          }
       }
    }
+   delete[] egval;
    egvec = NULL;
    printf("Done!\nNow to normalize the DOSs ..."); fflush(stdout);
  
@@ -1008,10 +1064,7 @@ void Phonon::ldos_egv()
  * ---------------------------------------------------------------------------- */
 void Phonon::ShowCell()
 {
-   printf("\n");
+   puts("==============================   Unit Cell Info   ==============================");
-   for (int i = 0; i < 30; ++i) printf("=");
-   printf("   Unit Cell Info   ");
-   for (int i = 0; i < 30; ++i) printf("="); printf("\n");
    printf("Number of atoms in the unit cell: %d\n", dynmat->nucell);
    printf("Basis  vectors  of the unit cell:\n");
    printf("  %15.8f  %15.8f  %15.8f\n", dynmat->basevec[0],  dynmat->basevec[1],  dynmat->basevec[2]);
@@ -1024,8 +1077,7 @@ void Phonon::ShowCell()
    printf("Atomic type and fractional coordinates:\n");
    for (int i = 0; i < dynmat->nucell; ++i)
      printf("%4d %12.8f %12.8f %12.8f\n", dynmat->attyp[i], dynmat->basis[i][0], dynmat->basis[i][1], dynmat->basis[i][2]);
-   for (int i = 0; i < 80; ++i) printf("=");
+   puts("================================================================================");
-   printf("\n");
  
    return;
 }
@@ -1101,7 +1153,7 @@ int Phonon::count_words(const char *line)
    strcpy(copy,line);
  
    char *ptr;
-   if (ptr = strchr(copy,'#')) *ptr = '\0';
+   if ((ptr = strchr(copy,'#'))) *ptr = '\0';
  
    if (strtok(copy," \t\n\r\f") == NULL) {
       memory->destroy(copy);

tools/phonon/phonon.h

@@ -1,22 +1,16 @@
 #ifndef PHONON_H
 #define PHONON_H
 
-#include "stdio.h"
-#include "stdlib.h"
-#include <complex>
-#include "dynmat.h"
-#include "memory.h"
-
-using namespace std;
-
 class Phonon{
 public:
-   Phonon(DynMat *);
+   Phonon(class DynMat *);
    ~Phonon();
  
-   DynMat *dynmat;
+   class DynMat *dynmat;
  
 private:
+   class UserInput *input;
+
    int nq, ndim, sysdim;
    double **qpts, *wt;
    double **eigs;
@@ -25,7 +19,7 @@ private:
    double *dos, fmin, fmax, df, rdf;
    double ***ldos;
  
-   Memory *memory;
+   class Memory *memory;
  
    void QMesh();
    void ComputeAll();

tools/phonon/phonopy.cpp

@@ -1,9 +1,20 @@
+
 #ifdef FFTW3
-#include <map>
+
 #include "phonopy.h"
-#include "math.h"
+
-#include "kpath.h"
+#include "global.h"
-#include "fftw3.h"
+#include "dynmat.h"
+#include "input.h"
+#include "memory.h"
+
+#include <fftw3.h>
+
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <map>
 
 /* ----------------------------------------------------------------------------
  * Class Phonopy is designed to interface with phonopy.
@@ -14,13 +25,16 @@ Phonopy::Phonopy(DynMat *dynmat)
    memory = new Memory();
    sysdim = dm->sysdim;
    fftdim = dm->fftdim;
+   input  = dm->input;
    fftdim2 = fftdim * fftdim;
    nucell = dm->nucell;
    nx = ny = nz = 5;
    write(1);
 
    char str[MAXLINE];
-   if (count_words(fgets(str,MAXLINE,stdin)) >= 3){
+   if (input == NULL) input = new UserInput(0);
+   input->read_stdin(str);
+   if (count_words(str) >= 3){
       nx = atoi(strtok(str," \t\n\r\f"));
       ny = atoi(strtok(NULL," \t\n\r\f"));
       nz = atoi(strtok(NULL," \t\n\r\f"));
@@ -36,7 +50,7 @@ Phonopy::Phonopy(DynMat *dynmat)
 
    memory->create(mass, nucell, "Phonopy:mass");
    for (int i = 0; i < nucell; ++i){
-      double m = 1./dm->M_inv_sqrt[i];
+      double m = 1.0/dm->M_inv_sqrt[i];
       mass[i] = m * m;
    }
 
@@ -68,7 +82,7 @@ return;
 void Phonopy::write(int flag)
 {
    if (flag == 1){  // basic information
-      for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+      puts("================================================================================");
       printf("Now to prepare the input files for phonopy.\n");
       printf("The dimension of your present supercell is   : %d x %d x %d.\n", dm->nx, dm->ny, dm->nz);
       printf("The size of the force constant matrix will be: %d x %d.\n", dm->npt*3, dm->npt*3);
@@ -84,19 +98,18 @@ void Phonopy::write(int flag)
    } else if (flag == 4){
       printf("Done!\nThe force constants information is extracted and written to FORCE_CONSTANTS,\n");
       printf("the primitive cell is written to POSCAR.primitive, and the input file for\n");
-      printf("phonopy band evaluation is written to band.conf.\n");
+      printf("phonopy band evaluation is written to band.conf.\n\n");
-      printf("One should be able to obtain the phonon band structure after correcting\n");
+      printf("One should be able to obtain the phonon band structure after\n");
-      printf("the element names in POSCAR.primitive and band.conf by running\n");
+      printf("  1) Correcting the `element names` in POSCAR.primitive and band.conf;\n");
-      printf("`phonopy --readfc -c POSCAR.primitive -p band.conf`.\n");
+      printf("  2) Running `phonopy --readfc -c POSCAR.primitive -p band.conf`.\n\n");
-      for (int ii = 0; ii < 80; ++ii) printf("-");
+      printf("Or the phonon density of states after\n");
-      printf("\n***          Remember to change the element names.           ***\n");
+      printf("  1) Correcting the `element names` in POSCAR.primitive and mesh.conf;\n");
-#ifdef UseSPG
+      printf("  2) Running `phonopy --readfc -c POSCAR.primitive -p mesh.conf`.\n");
-      for (int ii = 0; ii < 80; ++ii) printf("-");
+      puts("--------------------------------------------------------------------------------");
-#endif
+      printf("***         Remember to modify the `element names`.          ***\n");
 
    } else if (flag == 5){
-      for (int ii = 0; ii < 80; ++ii) printf("="); printf("\n");
+      puts("================================================================================");
-
    }
 return;
 }
@@ -162,7 +175,9 @@ void Phonopy::phonopy()
    memory->destroy(out);
 
    // in POSCAR, atoms are sorted/aggregated by type, while for LAMMPS there is no such requirment
-   int type_id[nucell], num_type[nucell], ntype = 0;
+   int *type_id = new int[nucell];
+   int *num_type = new int[nucell];
+   int ntype = 0;
    for (int i = 0; i < nucell; ++i) num_type[i] = 0;
    for (int i = 0; i < nucell; ++i){
       int ip = ntype;
@@ -221,14 +236,22 @@ void Phonopy::phonopy()
    // write the primitive cell in POSCAR format
    fp = fopen("POSCAR.primitive", "w");
    fprintf(fp, "Fix-phonon unit cell");
-   for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]);
+   for (int ip = 0; ip < ntype; ++ip){
+     for (int i = 0; i < nucell; ++i){
+       if (dm->attyp[i] == type_id[ip]){
+         fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]);
+         break;
+       }
+     }
+   }
    fprintf(fp, "\n1.\n"); 
    int ndim = 0;
    for (int idim = 0; idim < 3; ++idim){
       for (int jdim = 0; jdim < 3; ++jdim) fprintf(fp, "%lg ", dm->basevec[ndim++]);
       fprintf(fp, "\n");
    }
-   for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]); fprintf(fp, "\n");
+   for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]);
+   fprintf(fp, "\n");
    for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "%d ", num_type[ip]);
    fprintf(fp, "\nDirect\n");
    for (int ip = 0; ip < ntype; ++ip){
@@ -240,57 +263,48 @@ void Phonopy::phonopy()
    }
    fclose(fp);
 
-#ifdef UseSPG
+   // mesh.conf
-  // Get high symmetry k-path
+   fp = fopen("mesh.conf", "w");
-   QNodes *q = new QNodes();
+   fprintf(fp, "# From Fix-phonon");
-   kPath *kp = new kPath(dm, q);
+   for (int ip = 0; ip < ntype; ++ip){
-   kp->kpath();
+     for (int i = 0; i < nucell; ++i){
-#endif
+       if (dm->attyp[i] == type_id[ip]){
-   
+         fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]);
+         break;
+       }
+     }
+   }
+   fprintf(fp, "\n\nATOM_NAME = ");
+   for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]);
+   fprintf(fp, "\nDIM = %d %d %d\n", nx, ny, nz);
+   fprintf(fp, "MP  = 31 31 31\nFORCE_CONSTANTS = READ\n");
+   fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n");
+   fclose(fp);
+
+
    // band.conf
    fp = fopen("band.conf", "w");
    fprintf(fp, "# From Fix-phonon");
-   for (int ip = 0; ip < ntype; ++ip) fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[ip]);
+   for (int ip = 0; ip < ntype; ++ip){
+     for (int i = 0; i < nucell; ++i){
+       if (dm->attyp[i] == type_id[ip]){
+         fprintf(fp, ", Elem-%d: %lg", type_id[ip], mass[i]);
+         break;
+       }
+     }
+   }
    fprintf(fp, "\n\nATOM_NAME = ");
    for (int ip = 0; ip < ntype; ++ip) fprintf(fp, "Elem-%d ", type_id[ip]);
-   fprintf(fp, "\nDIM = %d %d %d\nBAND = ", nx, ny, nz);
+   fprintf(fp, "\nDIM = %d %d %d\nBAND = AUTO\n", nx, ny, nz);
-#ifdef UseSPG
+   fprintf(fp, "BAND_POINTS = 21\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n");
-   int nsect = q->qs.size();
+   fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n");
-   for (int i = 0; i < nsect; ++i){
-      fprintf(fp, " %lg %lg %lg", q->qs[i][0], q->qs[i][1], q->qs[i][2]);
-      if (i+1 < nsect){
-         double dq = 0.;
-         for (int j = 0; j < 3; ++j) dq += (q->qe[i][j] - q->qs[i+1][j]) * (q->qe[i][j] - q->qs[i+1][j]);
-         if (dq > ZERO) {
-            fprintf(fp, " %lg %lg %lg,", q->qe[i][0], q->qe[i][1], q->qe[i][2]);
-         }
-      } else if (i+1 == nsect){
-         fprintf(fp, " %lg %lg %lg\n", q->qe[i][0], q->qe[i][1], q->qe[i][2]);
-      }
-   }
-#endif
-   fprintf(fp, "\nBAND_POINTS = 21\nBAND_LABELS =");
-#ifdef UseSPG
-   for (int i = 0; i < q->ndstr.size(); ++i){
-      std::size_t found = q->ndstr[i].find("{/Symbol G}");
-      if (found != std::string::npos) q->ndstr[i].replace(found, found+11, "$\\Gamma$");
-      found = q->ndstr[i].find("/");
-      if (found != std::string::npos) q->ndstr[i].replace(found, found, " ");
-      fprintf(fp, " %s", q->ndstr[i].c_str());
-   }
-#endif
-   fprintf(fp, "\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n");
 
    // output info
    write(4);
-#ifdef UseSPG
-   kp->show_path();
-   delete kp;
-   delete q;
-#endif
    write(5);
-
+   delete[] type_id;
-return;
+   delete[] num_type;
+   return;
 }
 
 /*------------------------------------------------------------------------------
@@ -304,7 +318,7 @@ int Phonopy::count_words(const char *line)
    strcpy(copy,line);
   
    char *ptr;
-   if (ptr = strchr(copy,'#')) *ptr = '\0';
+   if ((ptr = strchr(copy,'#'))) *ptr = '\0';
   
    if (strtok(copy," \t\n\r\f") == NULL) {
      memory->destroy(copy);
@@ -314,7 +328,7 @@ int Phonopy::count_words(const char *line)
    while (strtok(NULL," \t\n\r\f")) n++;
   
    memory->destroy(copy);
-return n;
+   return n;
 }
 /*----------------------------------------------------------------------------*/
 #endif

tools/phonon/phonopy.h

@@ -4,22 +4,16 @@
 #ifndef PHONOPY_H
 #define PHONOPY_H
 
-#include "stdio.h"
+#include "zheevd.h"
-#include "stdlib.h"
-#include "string.h"
-#include "memory.h"
-#include "qnodes.h"
-#include "dynmat.h"
-#include "global.h"
 
 class Phonopy {
 public:
-   Phonopy(DynMat *);
+   Phonopy(class DynMat *);
    ~Phonopy();
 
 private:
-   Memory *memory;
+   class UserInput *input;
-   char str[MAXLINE];
+   class Memory *memory;
    int npt, fftdim2;       // local variables
    int nx, ny, nz, nucell; // local variables
    int sysdim, fftdim;     // local variables
@@ -27,7 +21,7 @@ private:
 
    doublecomplex **FC_all;
 
-   DynMat *dm;
+   class DynMat *dm;
    void write(int);
    void get_my_FC();
    void phonopy();

tools/phonon/timer.cpp

@@ -1,4 +1,7 @@
 #include "timer.h"
+
+#include <cstdio>
+
 /* -----------------------------------------------------------------------------
  * Initialization of time
  * -------------------------------------------------------------------------- */

tools/phonon/timer.h

@@ -1,9 +1,7 @@
 #ifndef TIMER_H
 #define TIMER_H
 
-#include "stdio.h"
+#include <ctime>
-#include "stdlib.h"
-#include "time.h"
 
 class Timer {
 public:

tools/phonon/tricubic/CMakeLists.txt

@@ -0,0 +1,18 @@
+
+# Support Linux from Ubuntu 20.04LTS onward, CentOS 7.x (with EPEL),
+# macOS, MSVC 2019 (=Version 16)
+cmake_minimum_required(VERSION 3.10)
+
+# set up project
+project(tricubic VERSION 1.1 DESCRIPTION "Tricubic library" LANGUAGES CXX)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+# hacks for MSVC to prevent lots of pointless warnings about "unsafe" functions
+if(MSVC)
+  add_compile_options(/wd4244)
+  add_compile_options(/wd4267)
+  add_compile_definitions(_CRT_SECURE_NO_WARNINGS)
+endif()
+
+add_library(tricubic STATIC tricubic.cpp)
+target_include_directories(tricubic PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

tools/phonon/tricubic/LICENSE

@@ -0,0 +1,340 @@
+		    GNU GENERAL PUBLIC LICENSE
+		       Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+     59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+			    Preamble
+
+  The licenses for most software are designed to take away your
+freedom to share and change it.  By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users.  This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it.  (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.)  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+  To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have.  You must make sure that they, too, receive or can get the
+source code.  And you must show them these terms so they know their
+rights.
+
+  We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+  Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software.  If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+  Finally, any free program is threatened constantly by software
+patents.  We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary.  To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+		    GNU GENERAL PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License.  The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language.  (Hereinafter, translation is included without limitation in
+the term "modification".)  Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope.  The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+  1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+  2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+    a) You must cause the modified files to carry prominent notices
+    stating that you changed the files and the date of any change.
+
+    b) You must cause any work that you distribute or publish, that in
+    whole or in part contains or is derived from the Program or any
+    part thereof, to be licensed as a whole at no charge to all third
+    parties under the terms of this License.
+
+    c) If the modified program normally reads commands interactively
+    when run, you must cause it, when started running for such
+    interactive use in the most ordinary way, to print or display an
+    announcement including an appropriate copyright notice and a
+    notice that there is no warranty (or else, saying that you provide
+    a warranty) and that users may redistribute the program under
+    these conditions, and telling the user how to view a copy of this
+    License.  (Exception: if the Program itself is interactive but
+    does not normally print such an announcement, your work based on
+    the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole.  If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works.  But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+  3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+    a) Accompany it with the complete corresponding machine-readable
+    source code, which must be distributed under the terms of Sections
+    1 and 2 above on a medium customarily used for software interchange; or,
+
+    b) Accompany it with a written offer, valid for at least three
+    years, to give any third party, for a charge no more than your
+    cost of physically performing source distribution, a complete
+    machine-readable copy of the corresponding source code, to be
+    distributed under the terms of Sections 1 and 2 above on a medium
+    customarily used for software interchange; or,
+
+    c) Accompany it with the information you received as to the offer
+    to distribute corresponding source code.  (This alternative is
+    allowed only for noncommercial distribution and only if you
+    received the program in object code or executable form with such
+    an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it.  For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable.  However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+  4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License.  Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+  5. You are not required to accept this License, since you have not
+signed it.  However, nothing else grants you permission to modify or
+distribute the Program or its derivative works.  These actions are
+prohibited by law if you do not accept this License.  Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+  6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions.  You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+  7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all.  For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices.  Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+  8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded.  In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+  9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number.  If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation.  If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+  10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission.  For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this.  Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+			    NO WARRANTY
+
+  11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+  12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+		     END OF TERMS AND CONDITIONS
+
+	    How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+    Gnomovision version 69, Copyright (C) year  name of author
+    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary.  Here is a sample; alter the names:
+
+  Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+  `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+  <signature of Ty Coon>, 1 April 1989
+  Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs.  If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library.  If this is what you want to do, use the GNU Library General
+Public License instead of this License.

tools/phonon/tricubic/README.md

@@ -0,0 +1,28 @@
+# libtricubic
+
+This folder contains a slightly refactored version of version 1.0 of
+the _libtricubic_ library, developed by François Lekien in 2004.
+
+The following paper explaining the method was published in 2005:
+
+> Lekien, F., & Marsden, J. (2005). _Tricubic interpolation in three
+> dimensions._ International Journal for Numerical Methods in Engineering,
+> 63(3), 455–471. doi:10.1002/nme.1296
+
+Some additional notes and the full matrix can be found in the technical notes:
+
+> Lekien, F., Coulliette, C., & Marsden, J. (2004). _Tricubic Engine - Technical
+> Notes and Full Matrix.
+
+The main article refers to the author's website
+(http://gyre.cds.caltech.edu/pub/software/tricubic/) to download the code and
+documentation. Unfortunately, this website no longer exists. Even the
+archive.org snapshot is useless; [A single snapshot from December 3rd 2009 is
+available](https://web.archive.org/web/20091203115835/http://gyre.cds.caltech.edu/pub/software/tricubic)
+which seems like an FTP listing. No files are accessible.
+
+The source code was obtained from https://github.com/nbigaouette/libtricubic/
+Only the sources for the library were retained. No functional changes were
+made, but some common programming conventions were applied, source files
+merged, and build support for CMake added.
+

tools/phonon/tricubic/tricubic.cpp

@@ -0,0 +1,185 @@
+
+#include "tricubic.h"
+
+#include <cmath>
+
+static const int A[64][64] = {
+  { 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 9,-9,-9, 9, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-6, 6, 6,-6, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-6, 6, 6,-6, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 4,-4,-4, 4, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3,-6,-3, 0, 0, 0, 0, 6,-6, 3,-3, 0, 0, 0, 0, 4, 2, 2, 1, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 3, 3, 0, 0, 0, 0,-4, 4,-2, 2, 0, 0, 0, 0,-2,-2,-1,-1, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 4, 2, 0, 0, 0, 0,-3, 3,-3, 3, 0, 0, 0, 0,-2,-1,-2,-1, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2,-2,-2, 0, 0, 0, 0, 2,-2, 2,-2, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0},
+  {-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 9,-9, 0, 0,-9, 9, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-6, 6, 0, 0, 6,-6, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0, 0, 0,-1, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 9,-9, 0, 0,-9, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 3, 0, 0,-6,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6,-6, 0, 0, 3,-3, 0, 0, 4, 2, 0, 0, 2, 1, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3,-3, 0, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 4, 0, 0,-2, 2, 0, 0,-2,-2, 0, 0,-1,-1, 0, 0},
+  { 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 9, 0,-9, 0,-9, 0, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 3, 0,-6, 0,-3, 0, 6, 0,-6, 0, 3, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 2, 0, 1, 0},
+  {-27,27,27,-27,27,-27,-27,27,-18,-9,18, 9,18, 9,-18,-9,-18,18,-9, 9,18,-18, 9,-9,-18,18,18,-18,-9, 9, 9,-9,-12,-6,-6,-3,12, 6, 6, 3,-12,-6,12, 6,-6,-3, 6, 3,-12,12,-6, 6,-6, 6,-3, 3,-8,-4,-4,-2,-4,-2,-2,-1},
+  {18,-18,-18,18,-18,18,18,-18, 9, 9,-9,-9,-9,-9, 9, 9,12,-12, 6,-6,-12,12,-6, 6,12,-12,-12,12, 6,-6,-6, 6, 6, 6, 3, 3,-6,-6,-3,-3, 6, 6,-6,-6, 3, 3,-3,-3, 8,-8, 4,-4, 4,-4, 2,-2, 4, 4, 2, 2, 2, 2, 1, 1},
+  {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 0,-3, 0, 3, 0, 3, 0,-4, 0, 4, 0,-2, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-2, 0,-1, 0,-1, 0},
+  {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6, 9,-9, 9,-9,-9, 9,-9, 9,12,-12,-12,12, 6,-6,-6, 6, 6, 3, 6, 3,-6,-3,-6,-3, 8, 4,-8,-4, 4, 2,-4,-2, 6,-6, 6,-6, 3,-3, 3,-3, 4, 2, 4, 2, 2, 1, 2, 1},
+  {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-6, 6,-6, 6, 6,-6, 6,-6,-8, 8, 8,-8,-4, 4, 4,-4,-3,-3,-3,-3, 3, 3, 3, 3,-4,-4, 4, 4,-2,-2, 2, 2,-4, 4,-4, 4,-2, 2,-2, 2,-2,-2,-2,-2,-1,-1,-1,-1},
+  { 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  {-6, 6, 0, 0, 6,-6, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 4,-4, 0, 0,-4, 4, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-6, 6, 0, 0, 6,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4,-2, 0, 0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-3, 3, 0, 0,-3, 3, 0, 0,-2,-1, 0, 0,-2,-1, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4,-4, 0, 0,-4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0,-2,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0},
+  {-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0,-6, 0, 6, 0, 6, 0,-6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-4, 0,-2, 0, 4, 0, 2, 0,-3, 0, 3, 0,-3, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0,-2, 0,-1, 0,-2, 0,-1, 0},
+  {18,-18,-18,18,-18,18,18,-18,12, 6,-12,-6,-12,-6,12, 6,12,-12, 6,-6,-12,12,-6, 6, 9,-9,-9, 9, 9,-9,-9, 9, 8, 4, 4, 2,-8,-4,-4,-2, 6, 3,-6,-3, 6, 3,-6,-3, 6,-6, 3,-3, 6,-6, 3,-3, 4, 2, 2, 1, 4, 2, 2, 1},
+  {-12,12,12,-12,12,-12,-12,12,-6,-6, 6, 6, 6, 6,-6,-6,-8, 8,-4, 4, 8,-8, 4,-4,-6, 6, 6,-6,-6, 6, 6,-6,-4,-4,-2,-2, 4, 4, 2, 2,-3,-3, 3, 3,-3,-3, 3, 3,-4, 4,-2, 2,-4, 4,-2, 2,-2,-2,-1,-1,-2,-2,-1,-1},
+  { 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0},
+  { 0, 0, 0, 0, 0, 0, 0, 0, 4, 0,-4, 0,-4, 0, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 2, 0,-2, 0,-2, 0, 2, 0,-2, 0, 2, 0,-2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0},
+  {-12,12,12,-12,12,-12,-12,12,-8,-4, 8, 4, 8, 4,-8,-4,-6, 6,-6, 6, 6,-6, 6,-6,-6, 6, 6,-6,-6, 6, 6,-6,-4,-2,-4,-2, 4, 2, 4, 2,-4,-2, 4, 2,-4,-2, 4, 2,-3, 3,-3, 3,-3, 3,-3, 3,-2,-1,-2,-1,-2,-1,-2,-1},
+  { 8,-8,-8, 8,-8, 8, 8,-8, 4, 4,-4,-4,-4,-4, 4, 4, 4,-4, 4,-4,-4, 4,-4, 4, 4,-4,-4, 4, 4,-4,-4, 4, 2, 2, 2, 2,-2,-2,-2,-2, 2, 2,-2,-2, 2, 2,-2,-2, 2,-2, 2,-2, 2,-2, 2,-2, 1, 1, 1, 1, 1, 1, 1, 1}
+};
+
+static const char tricubic_version_stored[] = "1.1";
+
+static int ijk2n(int i, int j, int k) {
+  return(i+4*j+16*k);
+}
+
+static void point2xyz(int p, int *x, int *y, int *z) {
+  switch (p) {
+  case 0: *x=0; *y=0; *z=0; break;
+  case 1: *x=1; *y=0; *z=0; break;
+  case 2: *x=0; *y=1; *z=0; break;
+  case 3: *x=1; *y=1; *z=0; break;
+  case 4: *x=0; *y=0; *z=1; break;
+  case 5: *x=1; *y=0; *z=1; break;
+  case 6: *x=0; *y=1; *z=1; break;
+  case 7: *x=1; *y=1; *z=1; break;
+  default:*x=0; *y=0; *z=0;
+  }
+}
+
+const char *tricubic_version(void) {
+  return(tricubic_version_stored);
+}
+
+void tricubic_pointID2xyz(int id, int *x, int *y, int *z) {
+  point2xyz(id,x,y,z);
+}
+
+void tricubic_pointID2xyz(int id, double *x, double *y, double *z) {
+  int x2,y2,z2;
+  point2xyz(id,&x2,&y2,&z2);
+  *x=(double)(x2);
+  *y=(double)(y2);
+  *z=(double)(z2);
+}
+
+void tricubic_get_coeff_stacked(double a[64], double x[64]) {
+  int i,j;
+  for (i=0;i<64;i++) {
+    a[i]=(double)(0.0);
+    for (j=0;j<64;j++) {
+      a[i]+=A[i][j]*x[j];
+    }
+  }
+}
+
+void tricubic_get_coeff(double a[64], double f[8], double dfdx[8], double dfdy[8], double dfdz[8], double d2fdxdy[8], double d2fdxdz[8], double d2fdydz[8], double d3fdxdydz[8]) {
+  int i;
+  double x[64];
+  for (i=0;i<8;i++) {
+    x[0+i]=f[i];
+    x[8+i]=dfdx[i];
+    x[16+i]=dfdy[i];
+    x[24+i]=dfdz[i];
+    x[32+i]=d2fdxdy[i];
+    x[40+i]=d2fdxdz[i];
+    x[48+i]=d2fdydz[i];
+    x[56+i]=d3fdxdydz[i];
+  }
+  tricubic_get_coeff_stacked(a,x);
+}
+
+double tricubic_eval(double a[64], double x, double y, double z) {
+  int i,j,k;
+  double ret=(double)(0.0);
+  /* TRICUBIC EVAL
+     This is the short version of tricubic_eval. It is used to compute
+     the value of the function at a given point (x,y,z). To compute
+     partial derivatives of f, use the full version with the extra args.
+  */
+  for (i=0;i<4;i++) {
+    for (j=0;j<4;j++) {
+      for (k=0;k<4;k++) {
+        ret+=a[ijk2n(i,j,k)]*pow(x,i)*pow(y,j)*pow(z,k);
+      }
+    }
+  }
+  return(ret);
+}
+
+double tricubic_eval(double a[64], double x, double y, double z, int derx, int dery, int derz) {
+  int i,j,k;
+  double ret=(double)(0.0);
+  double cont;
+  int w;
+  /* TRICUBIC_EVAL 
+     The full version takes 3 extra integers args that allows to evaluate
+     any partial derivative of f at the point
+     derx=dery=derz=0 => f
+     derx=2 dery=derz=0 => d2f/dx2
+     derx=dery=derz=1 =? d3f/dxdydz
+     NOTICE that (derx>3)||(dery>3)||(derz>3) => returns 0.0
+     this computes   \frac{\partial ^{derx+dery+derz} d}{\partial x ^{derx} \partial y ^{dery} \partial z ^{derz}}
+  */
+  for (i=derx;i<4;i++) {
+    for (j=dery;j<4;j++) {
+      for (k=derz;k<4;k++) {
+        cont=a[ijk2n(i,j,k)]*pow(x,i-derx)*pow(y,j-dery)*pow(z,k-derz);
+        for (w=0;w<derx;w++) {
+          cont*=(i-w);
+        }
+        for (w=0;w<dery;w++) {
+          cont*=(j-w);
+        }
+        for (w=0;w<derz;w++) {
+          cont*=(k-w);
+        }
+        ret+=cont;
+      }
+    }
+  }
+  return(ret);
+}

tools/phonon/tricubic/tricubic.h

@@ -0,0 +1,10 @@
+#ifndef TRICUBIC_H
+#define TRICUBIC_H
+extern const char *tricubic_version(void);
+extern void tricubic_get_coeff(double a[64], double f[8], double dfdx[8], double dfdy[8], double dfdz[8], double d2fdxdy[8], double d2fdxdz[8], double d2fdydz[8], double d3fdxdydz[8]);
+extern double tricubic_eval(double a[64], double x, double y, double z);
+extern double tricubic_eval(double a[64], double x, double y, double z, int derx, int dery, int derz);
+
+extern void tricubic_pointID2xyz(int id, int *x, int *y, int *z);
+extern void tricubic_pointID2xyz(int id, double *x, double *y, double *z);
+#endif

tools/phonon/version.h

@@ -1 +0,0 @@
-#define VERSION 21

tools/phonon/version.h.in

@@ -0,0 +1 @@
+#define VERSION @CMAKE_PROJECT_VERSION@

tools/phonon/zheevd.h

@@ -0,0 +1,16 @@
+#ifndef ZHEEVD_H
+#define ZHEEVD_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+  typedef struct { double r, i; } doublecomplex;
+
+  /* Subroutine */ int zheevd_(char *jobz, char *uplo, int *n,
+                               doublecomplex *a, int *lda, double *w, doublecomplex *work,
+                               int *lwork, double *rwork, int *lrwork, int *iwork,
+                               int *liwork, int *info);
+#ifdef __cplusplus
+}
+#endif
+#endif

unittest/commands/test_compute_global.cpp

@@ -294,6 +294,98 @@ TEST_F(ComputeGlobalTest, Reduction)
     EXPECT_DOUBLE_EQ(rep[2], 26);
     EXPECT_DOUBLE_EQ(rep[3], max[0]);
 }
+
+TEST_F(ComputeGlobalTest, Counts)
+{
+    if (lammps_get_natoms(lmp) == 0.0) GTEST_SKIP();
+
+    BEGIN_HIDE_OUTPUT();
+    command("pair_style zero 10.0");
+    command("pair_coeff * *");
+
+    command("variable t1 atom type==1");
+    command("variable t2 atom type==2");
+    command("variable t3 atom type==3");
+    command("variable t4 atom type==4");
+    command("variable t5 atom type==5");
+    command("compute tsum all reduce sum v_t1 v_t2 v_t3 v_t4 v_t5");
+    command("compute tcnt all count/type atom");
+    command("compute bcnt all count/type bond");
+    command("compute acnt all count/type angle");
+    command("compute dcnt all count/type dihedral");
+    command("compute icnt all count/type improper");
+    command("thermo_style custom c_tsum[*] c_tcnt[*] c_bcnt[*] c_acnt[*] c_dcnt[*] c_icnt[*]");
+    command("run 0 post no");
+    END_HIDE_OUTPUT();
+
+    auto tsum = get_vector("tsum");
+    auto tcnt = get_vector("tcnt");
+    auto bcnt = get_vector("bcnt");
+    auto bbrk = get_scalar("bcnt");
+    auto acnt = get_vector("acnt");
+    auto dcnt = get_vector("dcnt");
+    auto icnt = get_vector("icnt");
+
+    EXPECT_DOUBLE_EQ(tsum[0], tcnt[0]);
+    EXPECT_DOUBLE_EQ(tsum[1], tcnt[1]);
+    EXPECT_DOUBLE_EQ(tsum[2], tcnt[2]);
+    EXPECT_DOUBLE_EQ(tsum[3], tcnt[3]);
+    EXPECT_DOUBLE_EQ(tsum[4], tcnt[4]);
+
+    EXPECT_DOUBLE_EQ(bbrk, 0.0);
+
+    EXPECT_DOUBLE_EQ(bcnt[0], 3.0);
+    EXPECT_DOUBLE_EQ(bcnt[1], 6.0);
+    EXPECT_DOUBLE_EQ(bcnt[2], 3.0);
+    EXPECT_DOUBLE_EQ(bcnt[3], 2.0);
+    EXPECT_DOUBLE_EQ(bcnt[4], 10.0);
+
+    EXPECT_DOUBLE_EQ(acnt[0], 6.0);
+    EXPECT_DOUBLE_EQ(acnt[1], 10.0);
+    EXPECT_DOUBLE_EQ(acnt[2], 5.0);
+    EXPECT_DOUBLE_EQ(acnt[3], 9.0);
+
+    EXPECT_DOUBLE_EQ(dcnt[0], 3.0);
+    EXPECT_DOUBLE_EQ(dcnt[1], 8.0);
+    EXPECT_DOUBLE_EQ(dcnt[2], 3.0);
+    EXPECT_DOUBLE_EQ(dcnt[3], 4.0);
+    EXPECT_DOUBLE_EQ(dcnt[4], 13.0);
+
+    EXPECT_DOUBLE_EQ(icnt[0], 1.0);
+    EXPECT_DOUBLE_EQ(icnt[1], 1.0);
+
+    BEGIN_HIDE_OUTPUT();
+    command("delete_bonds all bond 3 remove");
+    command("run 0 post no");
+    END_HIDE_OUTPUT();
+
+    bcnt = get_vector("bcnt");
+    bbrk = get_scalar("bcnt");
+    acnt = get_vector("acnt");
+    dcnt = get_vector("dcnt");
+    icnt = get_vector("icnt");
+
+    EXPECT_DOUBLE_EQ(bbrk, 0.0);
+    EXPECT_DOUBLE_EQ(bcnt[0], 3.0);
+    EXPECT_DOUBLE_EQ(bcnt[1], 6.0);
+    EXPECT_DOUBLE_EQ(bcnt[2], 0.0);
+    EXPECT_DOUBLE_EQ(bcnt[3], 2.0);
+    EXPECT_DOUBLE_EQ(bcnt[4], 10.0);
+
+    EXPECT_DOUBLE_EQ(acnt[0], 6.0);
+    EXPECT_DOUBLE_EQ(acnt[1], 10.0);
+    EXPECT_DOUBLE_EQ(acnt[2], 5.0);
+    EXPECT_DOUBLE_EQ(acnt[3], 9.0);
+
+    EXPECT_DOUBLE_EQ(dcnt[0], 3.0);
+    EXPECT_DOUBLE_EQ(dcnt[1], 8.0);
+    EXPECT_DOUBLE_EQ(dcnt[2], 3.0);
+    EXPECT_DOUBLE_EQ(dcnt[3], 4.0);
+    EXPECT_DOUBLE_EQ(dcnt[4], 13.0);
+
+    EXPECT_DOUBLE_EQ(icnt[0], 1.0);
+    EXPECT_DOUBLE_EQ(icnt[1], 1.0);
+}
 } // namespace LAMMPS_NS
 
 int main(int argc, char **argv)