diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index a71347c2c4..88b7ec422e 100644 --- a/cmake/CMakeLists.txt +++ b/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) diff --git a/cmake/Modules/Tools.cmake b/cmake/Modules/Tools.cmake index c4c33cc40c..285c5f2405 100644 --- a/cmake/Modules/Tools.cmake +++ b/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) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 2f9c8b03fb..755000c976 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -46,6 +46,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`com/chunk ` * :doc:`contact/atom ` * :doc:`coord/atom (k) ` + * :doc:`count/type ` * :doc:`damage/atom ` * :doc:`dihedral ` * :doc:`dihedral/local ` diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index 88c8acb87f..9fcd676c36 100644 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -1,48 +1,56 @@ Broken Bonds ============ -Typically, bond interactions persist for the duration of a simulation in -LAMMPS. However, there are some exceptions that allow for bonds to -break, including the :doc:`quartic bond style ` and the -bond styles in the :doc:`BPM package ` which contains the -:doc:`bpm/spring ` and :doc:`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. +Typically, molecular bond interactions persist for the duration of a +simulation in LAMMPS. However, some commands break bonds dynamically, +including the following: -Users are normally able to weight the contribution of pair forces to atoms -that are bonded using the :doc:`special_bonds command `. -When bonds break, this is not always the case. For the quartic bond style, -pair forces are always turned off between bonded particles. LAMMPS does -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 ` command, whether the bond -is broken or not. This means that :doc:`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. +* :doc:`bond_style quartic ` +* :doc:`fix bond/break ` +* :doc:`fix bond/react ` +* :doc:`BPM package ` bond styles -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 ` -page. +A bond can break if it is stretched beyond a user-defined threshold or +more generally if other criteria are met. -Bonds can also be broken by fixes which change bond topology, including -:doc:`fix bond/break ` and -:doc:`fix bond/react `. These fixes will automatically -trigger a rebuild of the neighbor list and update special bond data structures -when bonds are broken. +For the quartic bond style, when a bond is broken its bond type is set +to 0 to effectively break it and pairwise forces between the two atoms +in the broken bond are "turned on". Angles, dihedrals, etc cannot be +defined for a system when :doc:`bond_style quartic ` is +used. -Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. The -:doc:`delete_bonds ` command can also be used to query the -status of broken bonds or permanently delete them, e.g.: +Similarly, bond styles in the BPM package are also incompatible with +angles, dihedrals, etc. and when a bond breaks its type is set to zero. +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 +` page. + +The :doc:`fix bond/break ` and :doc:`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 +` command, bonds with type 0 are not included. + +The :doc:`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 ` can also be used -to tally the current number of bonds per atom, excluding broken bonds. +The compute :doc:`count/type ` command tallies the +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 ` command tallies the +current number of bonds each atom is part of, excluding broken bonds +with type = 0. diff --git a/doc/src/Tools.rst b/doc/src/Tools.rst index 2eae1aabc7..bccb9abe73 100644 --- a/doc/src/Tools.rst +++ b/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 -analyzing the output of the :doc:`fix phonon ` command in -the PHONON package. +The phonon subdirectory contains a post-processing tool, *phana*, useful +for analyzing the output of the :doc:`fix phonon ` command +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 diff --git a/doc/src/bond_style.rst b/doc/src/bond_style.rst index b33d0a9e9a..95f463e695 100644 --- a/doc/src/bond_style.rst +++ b/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 ` command. Bonds are defined between specified pairs of atoms and -remain in force for the duration of the simulation (unless the bond -breaks which is possible in some bond potentials). The list of bonded -atoms is read in by a :doc:`read_data ` or -:doc:`read_restart ` command from a data or restart file. -By contrast, pair potentials are typically defined between all pairs -of atoms within a cutoff distance and the set of active interactions -changes over time. +remain in force for the duration of the simulation (unless new bonds +are created or existing bonds break, which is possible in some fixes +and bond potentials). The list of bonded atoms is read in by a +:doc:`read_data ` or :doc:`read_restart ` +command from a data or restart file. By contrast, pair potentials are +typically defined between all pairs of atoms within a cutoff distance +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. diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 880f60a8a6..950487cb72 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -200,6 +200,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`com/chunk ` - center of mass for each chunk * :doc:`contact/atom ` - contact count for each spherical particle * :doc:`coord/atom ` - coordination number for each atom +* :doc:`count/type ` - count of atoms or bonds by type * :doc:`damage/atom ` - Peridynamic damage for each atom * :doc:`dihedral ` - energy of each dihedral sub-style * :doc:`dihedral/local ` - angle of each dihedral diff --git a/doc/src/compute_count_type.rst b/doc/src/compute_count_type.rst new file mode 100644 index 0000000000..73e9d50f73 --- /dev/null +++ b/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 ` 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 `, +:doc:`fix pour `, or :doc:`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 +` command or defined by the :doc:`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 ` +* :doc:`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 +` and :doc:`BPM bond styles ` are exceptions, +see the discussion below). Thus they are not included in the counts +for each type: + +* :doc:`delete_bonds remove ` +* :doc:`bond_style quartic ` +* :doc:`fix bond/react ` +* :doc:`fix bond/create ` +* :doc:`fix bond/break ` +* :doc:`BPM package ` 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 ` and :doc:`BPM +bond styles ` break bonds by doing this. See the :doc:` +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 +` 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 diff --git a/src/KOKKOS/fix_langevin_kokkos.cpp b/src/KOKKOS/fix_langevin_kokkos.cpp index 96fa58f601..4d7a3e8820 100644 --- a/src/KOKKOS/fix_langevin_kokkos.cpp +++ b/src/KOKKOS/fix_langevin_kokkos.cpp @@ -44,6 +44,7 @@ FixLangevinKokkos::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::initial_integrate_item(int i) const /* ---------------------------------------------------------------------- */ +template +void FixLangevinKokkos::fused_integrate(int vflag) +{ + initial_integrate(vflag); +} + +/* ---------------------------------------------------------------------- */ + template void FixLangevinKokkos::post_force(int /*vflag*/) { diff --git a/src/KOKKOS/fix_langevin_kokkos.h b/src/KOKKOS/fix_langevin_kokkos.h index 3f23f57cc3..4fc22a1df1 100644 --- a/src/KOKKOS/fix_langevin_kokkos.h +++ b/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; diff --git a/src/KOKKOS/fix_nve_kokkos.cpp b/src/KOKKOS/fix_nve_kokkos.cpp index 5dcb611d41..59cc90c088 100644 --- a/src/KOKKOS/fix_nve_kokkos.cpp +++ b/src/KOKKOS/fix_nve_kokkos.cpp @@ -29,6 +29,7 @@ FixNVEKokkos::FixNVEKokkos(LAMMPS *lmp, int narg, char **arg) : FixNVE(lmp, narg, arg) { kokkosable = 1; + fuse_integrate_flag = 1; atomKK = (AtomKokkos *) atom; execution_space = ExecutionSpaceFromDevice::space; @@ -159,6 +160,66 @@ void FixNVEKokkos::final_integrate_rmass_item(int i) const } } +/* ---------------------------------------------------------------------- + allow for both per-type and per-atom mass +------------------------------------------------------------------------- */ + +template +void FixNVEKokkos::fused_integrate(int /*vflag*/) +{ + atomKK->sync(execution_space,datamask_read); + + x = atomKK->k_x.view(); + v = atomKK->k_v.view(); + f = atomKK->k_f.view(); + rmass = atomKK->k_rmass.view(); + mass = atomKK->k_mass.view(); + type = atomKK->k_type.view(); + mask = atomKK->k_mask.view(); + int nlocal = atomKK->nlocal; + if (igroup == atomKK->firstgroup) nlocal = atomKK->nfirst; + + if (rmass.data()) { + FixNVEKokkosFusedIntegrateFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } else { + FixNVEKokkosFusedIntegrateFunctor functor(this); + Kokkos::parallel_for(nlocal,functor); + } + + atomKK->modified(execution_space,datamask_modify); +} + +template +KOKKOS_INLINE_FUNCTION +void FixNVEKokkos::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 +KOKKOS_INLINE_FUNCTION +void FixNVEKokkos::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 @@ -168,6 +229,8 @@ void FixNVEKokkos::cleanup_copy() vatom = nullptr; } +/* ---------------------------------------------------------------------- */ + namespace LAMMPS_NS { template class FixNVEKokkos; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_kokkos.h b/src/KOKKOS/fix_nve_kokkos.h index 215aacb4a0..a1e8e1398c 100644 --- a/src/KOKKOS/fix_nve_kokkos.h +++ b/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 +struct FixNVEKokkosFusedIntegrateFunctor { + typedef DeviceType device_type ; + FixNVEKokkos c; + + FixNVEKokkosFusedIntegrateFunctor(FixNVEKokkos* 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 diff --git a/src/KOKKOS/fix_nve_sphere_kokkos.cpp b/src/KOKKOS/fix_nve_sphere_kokkos.cpp index 51e57b839e..38f6a40792 100644 --- a/src/KOKKOS/fix_nve_sphere_kokkos.cpp +++ b/src/KOKKOS/fix_nve_sphere_kokkos.cpp @@ -28,6 +28,7 @@ FixNVESphereKokkos::FixNVESphereKokkos(LAMMPS *lmp, int narg, char * FixNVESphere(lmp, narg, arg) { kokkosable = 1; + fuse_integrate_flag = 1; atomKK = (AtomKokkos *)atom; execution_space = ExecutionSpaceFromDevice::space; @@ -164,6 +165,73 @@ void FixNVESphereKokkos::final_integrate_item(const int i) const } } +/* ---------------------------------------------------------------------- */ + +template +void FixNVESphereKokkos::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(); + v = atomKK->k_v.view(); + omega = atomKK->k_omega.view(); + f = atomKK->k_f.view(); + torque = atomKK->k_torque.view(); + mask = atomKK->k_mask.view(); + rmass = atomKK->k_rmass.view(); + radius = atomKK->k_radius.view(); + mu = atomKK->k_mu.view(); + + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + FixNVESphereKokkosFusedIntegrateFunctor 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 +KOKKOS_INLINE_FUNCTION +void FixNVESphereKokkos::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; #ifdef LMP_KOKKOS_GPU diff --git a/src/KOKKOS/fix_nve_sphere_kokkos.h b/src/KOKKOS/fix_nve_sphere_kokkos.h index 268b4ea9ce..02acb466f2 100644 --- a/src/KOKKOS/fix_nve_sphere_kokkos.h +++ b/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::t_x_array x; @@ -77,6 +80,17 @@ struct FixNVESphereKokkosFinalIntegrateFunctor { } }; +template +struct FixNVESphereKokkosFusedIntegrateFunctor { + typedef DeviceType device_type; + FixNVESphereKokkos c; + FixNVESphereKokkosFusedIntegrateFunctor(FixNVESphereKokkos *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 diff --git a/src/KOKKOS/modify_kokkos.cpp b/src/KOKKOS/modify_kokkos.cpp index 9d8c16603e..0b81a1cabb 100644 --- a/src/KOKKOS/modify_kokkos.cpp +++ b/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; +} diff --git a/src/KOKKOS/modify_kokkos.h b/src/KOKKOS/modify_kokkos.h index 5edf5cd662..527518219c 100644 --- a/src/KOKKOS/modify_kokkos.h +++ b/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: }; diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h index 0ff244f67d..c244d6a8f9 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/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; diff --git a/src/KOKKOS/verlet_kokkos.cpp b/src/KOKKOS/verlet_kokkos.cpp index 01401542cd..7570f1d8fa 100644 --- a/src/KOKKOS/verlet_kokkos.cpp +++ b/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; +} diff --git a/src/KOKKOS/verlet_kokkos.h b/src/KOKKOS/verlet_kokkos.h index 067df54f4f..c71211c542 100644 --- a/src/KOKKOS/verlet_kokkos.h +++ b/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); }; } diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp new file mode 100644 index 0000000000..3cb2e6bc3a --- /dev/null +++ b/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; +} diff --git a/src/compute_count_type.h b/src/compute_count_type.h new file mode 100644 index 0000000000..b766edfae8 --- /dev/null +++ b/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 diff --git a/src/fix.cpp b/src/fix.cpp index 78fdf1ce65..82389ba433 100644 --- a/src/fix.cpp +++ b/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; } /* ---------------------------------------------------------------------- */ diff --git a/src/fix.h b/src/fix.h index 334f61ff2b..30373ab6f2 100644 --- a/src/fix.h +++ b/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 *) {} diff --git a/src/lmptype.h b/src/lmptype.h index bf56d07d77..ecd6ee5761 100644 --- a/src/lmptype.h +++ b/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 diff --git a/src/modify.h b/src/modify.h index 7a3f54c277..6ca4b4ad26 100644 --- a/src/modify.h +++ b/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(); diff --git a/src/pair.cpp b/src/pair.cpp index 9ae24e1e93..34c8e4978e 100644 --- a/src/pair.cpp +++ b/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; - reverse_comm_device = 0; - copymode = 0; + kokkosable = copymode = 0; + reverse_comm_device = fuse_force_clear_flag = 0; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair.h b/src/pair.h index 8c856660e9..d36b48f340 100644 --- a/src/pair.h +++ b/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; diff --git a/src/timer.h b/src/timer.h index 5c100db1c0..f7efa5ac64 100644 --- a/src/timer.h +++ b/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); } diff --git a/tools/phonon/CMakeLists.spglib b/tools/phonon/CMakeLists.spglib new file mode 100644 index 0000000000..566e58a2b6 --- /dev/null +++ b/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) + diff --git a/tools/phonon/CMakeLists.txt b/tools/phonon/CMakeLists.txt new file mode 100644 index 0000000000..60da1cc79f --- /dev/null +++ b/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 $) + +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 "$") + set(LAPACK_LIBRARIES "$") + 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}) diff --git a/tools/phonon/Makefile b/tools/phonon/Makefile deleted file mode 100644 index 6afae5f312..0000000000 --- a/tools/phonon/Makefile +++ /dev/null @@ -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 $< diff --git a/tools/phonon/README b/tools/phonon/README index db69ac50c8..c1430f5940 100644 --- a/tools/phonon/README +++ b/tools/phonon/README @@ -5,34 +5,38 @@ analyse the phonon related information. #------------------------------------------------------------------------------- 1. Dependencies - The clapack library is needed to solve the eigen problems, - which could be downloaded from: - http://www.netlib.org/clapack/ - - The tricubic library is also needed to do tricubic interpolations, - which could now be obtained from: - https://github.com/nbigaouette/libtricubic/ - + The ZHEEVD LAPACK function is needed to solve the eigen problems. + A C++ compilable version based on CLAPACK is included in the linalg folder + and will be automatically built. + + The tricubic library is also needed to do tricubic interpolations. + A copy is included and will be automatically built. + 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. - It can be obtained from: - http://spglib.sourceforge.net/ + automatic mode. Currently, version 1.11.2.1 of spglib is used. + It is automatically downloaded and compiled unless the -DUSE_SPGLIB=off + 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: - http://www.fftw.org - + It is autodetected and used if available. + + FFTW 3 can be downloaded from: http://www.fftw.org + 2. Compilation - To compile the code, one needs therefore to install the above - libraries and set the paths correctly in the Makefile. - Once this is done, by typing - make - will yield the executable "phana". + To compile the code, one needs to have CMake version 3.16 + or later installed. + + The CMake configuration is done with: + 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 diff --git a/tools/phonon/disp.cpp b/tools/phonon/disp.cpp index 8a53873383..79a22aeee4 100644 --- a/tools/phonon/disp.cpp +++ b/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 +#include +#include +#include +#include +#include /*------------------------------------------------------------------------------ * 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 (count_words(fgets(str, MAXLINE, stdin)) < 3); + while ( 1 ){ + 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; diff --git a/tools/phonon/dynmat.cpp b/tools/phonon/dynmat.cpp index a6c4105547..aeb8fe6430 100644 --- a/tools/phonon/dynmat.cpp +++ b/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 +#include +#include /* ---------------------------------------------------------------------------- * 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"); - for (int i = 0; i < 80; ++i) printf("="); printf("\n\n"); - + delete[] egvs; + 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; } /* --------------------------------------------------------------------*/ diff --git a/tools/phonon/dynmat.h b/tools/phonon/dynmat.h index 10caa1ce2a..e11badfde6 100644 --- a/tools/phonon/dynmat.h +++ b/tools/phonon/dynmat.h @@ -1,11 +1,9 @@ #ifndef DYNMAT_H #define DYNMAT_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "interpolate.h" +#include "zheevd.h" + +#include 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; - - Memory *memory; + class Interpolate *interpolate; + class Memory *memory; - int nasr; void EnforceASR(); char *binfile, *dmfile; - double boltz, q[3]; + double boltz; doublecomplex **DM_all; diff --git a/tools/phonon/global.h b/tools/phonon/global.h index c20a91f317..bd1608ee07 100644 --- a/tools/phonon/global.h +++ b/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)) diff --git a/tools/phonon/green.cpp b/tools/phonon/green.cpp index 41f5b14886..0f687fd834 100644 --- a/tools/phonon/green.cpp +++ b/tools/phonon/green.cpp @@ -1,10 +1,10 @@ -#include -#include -#include -#include #include "green.h" + +#include "memory.h" + #include -#include "global.h" +#include +#include /******************************************************************************* * 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 Z, rec_x, rec_x_inv; - std::complex cunit = std::complex(0.,1.); - double w = wmin; for (int i = 0; i < nw; ++i){ diff --git a/tools/phonon/green.h b/tools/phonon/green.h index 1a137e53ba..21e3b091dc 100644 --- a/tools/phonon/green.h +++ b/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 diff --git a/tools/phonon/input.cpp b/tools/phonon/input.cpp new file mode 100644 index 0000000000..c2059043c7 --- /dev/null +++ b/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; +} +/* ---------------------------------------------------------------- */ diff --git a/tools/phonon/input.h b/tools/phonon/input.h new file mode 100644 index 0000000000..c931895c4d --- /dev/null +++ b/tools/phonon/input.h @@ -0,0 +1,17 @@ +#ifndef INPUT_H +#define INPUT_H + +#include + +class UserInput { +public: + UserInput(int); + ~UserInput(); + + void read_stdin(char *); + +private: + FILE *fp; + +}; +#endif diff --git a/tools/phonon/interpolate.cpp b/tools/phonon/interpolate.cpp index 09e261c763..8ea551d1a8 100644 --- a/tools/phonon/interpolate.cpp +++ b/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 +#include +#include /* ---------------------------------------------------------------------------- * 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("="); - printf("\nWhich interpolation method would you like to use?\n"); + if (input == NULL) input = new UserInput(0); + + 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(); diff --git a/tools/phonon/interpolate.h b/tools/phonon/interpolate.h index c650b30908..b9e0242b96 100644 --- a/tools/phonon/interpolate.h +++ b/tools/phonon/interpolate.h @@ -1,16 +1,7 @@ #ifndef INTERPOLATION_H #define INTERPOLATION_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "tricubic.h" - -extern "C"{ -#include "f2c.h" -#include "clapack.h" -} +#include "zheevd.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; diff --git a/tools/phonon/kpath.cpp b/tools/phonon/kpath.cpp index 842e680611..49730b42b6 100644 --- a/tools/phonon/kpath.cpp +++ b/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 +#include +#include +#include /* ---------------------------------------------------------------------------- * 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; diff --git a/tools/phonon/kpath.h b/tools/phonon/kpath.h index bbcedf15d4..ade1e3a27c 100644 --- a/tools/phonon/kpath.h +++ b/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(DynMat *, QNodes *); + kPath(class DynMat *, class QNodes *); ~kPath(); void kpath(); @@ -19,13 +14,11 @@ public: void show_info(); private: - - Memory *memory; - - DynMat *dynmat; - QNodes *q; + class Memory *memory; + class DynMat *dynmat; + class QNodes *q; char symbol[11]; - int spgnum, sysdim, fftdim, num_atom, *attyp; + int spgnum, sysdim, num_atom, *attyp; double latvec[3][3], **atpos; }; diff --git a/tools/phonon/main.cpp b/tools/phonon/main.cpp index d7d5baa1cc..d0193a037c 100644 --- a/tools/phonon/main.cpp +++ b/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) { diff --git a/tools/phonon/memory.cpp b/tools/phonon/memory.cpp index 4d65e83ef8..18dfa2fa43 100644 --- a/tools/phonon/memory.cpp +++ b/tools/phonon/memory.cpp @@ -1,8 +1,8 @@ -#include "stdio.h" -#include "stdlib.h" -#include "string.h" #include "memory.h" +#include +#include + /* ---------------------------------------------------------------------- safe malloc ------------------------------------------------------------------------- */ diff --git a/tools/phonon/memory.h b/tools/phonon/memory.h index ae2feceba3..13eeca4b14 100644 --- a/tools/phonon/memory.h +++ b/tools/phonon/memory.h @@ -4,18 +4,16 @@ #define __STDC_LIMIT_MACROS #define __STDC_FORMAT_MACROS -#include "stdio.h" -#include "stdlib.h" -#include "limits.h" -#include "stdint.h" -#include "inttypes.h" +#include +#include +#include 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - 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 - 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]; - n += n2; + array[i] = &data[n]; + n += n2; } return array; } template - 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 - 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 - cannot grow it -------------------------------------------------------------------------- */ + create a 3d array + ------------------------------------------------------------------------- */ template - TYPE **create2d_offset(TYPE **&array, int n1, int n2lo, int n2hi, - 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 - 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 - void destroy2d_offset(TYPE **array, int offset) - { - if (array == NULL) return; - sfree(&array[0][offset]); - sfree(array); - } - -/* ---------------------------------------------------------------------- - create a 3d array -------------------------------------------------------------------------- */ - - template - TYPE ***create(TYPE ***&array, int n1, int n2, int n3, const char *name) + 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]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - 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 - 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]; - for (j = 0; j < n2; j++) { - plane[i*n2+j] = &data[n]; - n += n3; - } + array[i] = &plane[i*n2]; + for (j = 0; j < n2; j++) { + plane[i*n2+j] = &data[n]; + n += n3; + } } return array; } template - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - TYPE *****create(TYPE *****&array, int n1, int n2, int n3, int n4, - const char *name) - {fail(name);} - -/* ---------------------------------------------------------------------- - destroy a 4d array -------------------------------------------------------------------------- */ - - template - 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 - bigint usage(TYPE *array, int n) + bigint usage(TYPE *, int n) { bigint bytes = sizeof(TYPE) * n; return bytes; } template - 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 - 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 - 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 diff --git a/tools/phonon/phonon.cpp b/tools/phonon/phonon.cpp index 56fa409e06..06372dcd1b 100644 --- a/tools/phonon/phonon.cpp +++ b/tools/phonon/phonon.cpp @@ -1,9 +1,18 @@ -#include -#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 +#include +#include +#include #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"); - for (int i = 0; i < 37; ++i) printf("="); - printf(" Menu "); - for (int i = 0; i < 37; ++i) printf("="); printf("\n"); + puts("\n===================================== Menu ====================================="); 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 (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + 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 (count_words(fgets(str,MAXLINE,stdin)) < 3); + while ( 1 ){ + 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]; - double symprec = 1.e-4, pos[num_atom][3]; + int **grid_point; + 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; - - double wsum = 0.; + delete[] iq2idx; + delete[] map; + 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"); - for (int i = 0; i < 30; ++i) printf("="); - printf(" Unit Cell Info "); - for (int i = 0; i < 30; ++i) printf("="); printf("\n"); + puts("============================== Unit Cell Info =============================="); 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("="); - printf("\n"); + puts("================================================================================"); 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); diff --git a/tools/phonon/phonon.h b/tools/phonon/phonon.h index 69a1fe5d50..46ce17f268 100644 --- a/tools/phonon/phonon.h +++ b/tools/phonon/phonon.h @@ -1,22 +1,16 @@ #ifndef PHONON_H #define PHONON_H -#include "stdio.h" -#include "stdlib.h" -#include -#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(); diff --git a/tools/phonon/phonopy.cpp b/tools/phonon/phonopy.cpp index cfe32ab61a..2cee319aa7 100644 --- a/tools/phonon/phonopy.cpp +++ b/tools/phonon/phonopy.cpp @@ -1,9 +1,20 @@ + #ifdef FFTW3 -#include + #include "phonopy.h" -#include "math.h" -#include "kpath.h" -#include "fftw3.h" + +#include "global.h" +#include "dynmat.h" +#include "input.h" +#include "memory.h" + +#include + +#include +#include +#include +#include +#include /* ---------------------------------------------------------------------------- * 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("One should be able to obtain the phonon band structure after correcting\n"); - printf("the element names in POSCAR.primitive and band.conf by running\n"); - printf("`phonopy --readfc -c POSCAR.primitive -p band.conf`.\n"); - for (int ii = 0; ii < 80; ++ii) printf("-"); - printf("\n*** Remember to change the element names. ***\n"); -#ifdef UseSPG - for (int ii = 0; ii < 80; ++ii) printf("-"); -#endif + printf("phonopy band evaluation is written to band.conf.\n\n"); + printf("One should be able to obtain the phonon band structure after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and band.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p band.conf`.\n\n"); + printf("Or the phonon density of states after\n"); + printf(" 1) Correcting the `element names` in POSCAR.primitive and mesh.conf;\n"); + printf(" 2) Running `phonopy --readfc -c POSCAR.primitive -p mesh.conf`.\n"); + puts("--------------------------------------------------------------------------------"); + 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 - // Get high symmetry k-path - QNodes *q = new QNodes(); - kPath *kp = new kPath(dm, q); - kp->kpath(); -#endif - + // mesh.conf + fp = fopen("mesh.conf", "w"); + fprintf(fp, "# From Fix-phonon"); + 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\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); -#ifdef UseSPG - int nsect = q->qs.size(); - 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"); + fprintf(fp, "\nDIM = %d %d %d\nBAND = AUTO\n", nx, ny, nz); + fprintf(fp, "BAND_POINTS = 21\nFORCE_CONSTANTS = READ\nBAND_CONNECTION = .TRUE.\n"); + fprintf(fp, "#FC_SYMMETRY = .TRUE.\n#SYMMETRY_TOLERANCE = 0.01\n"); // output info write(4); -#ifdef UseSPG - kp->show_path(); - delete kp; - delete q; -#endif write(5); - -return; + delete[] type_id; + 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 diff --git a/tools/phonon/phonopy.h b/tools/phonon/phonopy.h index 2f3006e791..e1200b8242 100644 --- a/tools/phonon/phonopy.h +++ b/tools/phonon/phonopy.h @@ -4,22 +4,16 @@ #ifndef PHONOPY_H #define PHONOPY_H -#include "stdio.h" -#include "stdlib.h" -#include "string.h" -#include "memory.h" -#include "qnodes.h" -#include "dynmat.h" -#include "global.h" +#include "zheevd.h" class Phonopy { public: - Phonopy(DynMat *); + Phonopy(class DynMat *); ~Phonopy(); private: - Memory *memory; - char str[MAXLINE]; + class UserInput *input; + 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(); diff --git a/tools/phonon/timer.cpp b/tools/phonon/timer.cpp index 953257340d..b94bed40b1 100644 --- a/tools/phonon/timer.cpp +++ b/tools/phonon/timer.cpp @@ -1,4 +1,7 @@ #include "timer.h" + +#include + /* ----------------------------------------------------------------------------- * Initialization of time * -------------------------------------------------------------------------- */ diff --git a/tools/phonon/timer.h b/tools/phonon/timer.h index cd04dea56d..993a33f9f3 100644 --- a/tools/phonon/timer.h +++ b/tools/phonon/timer.h @@ -1,9 +1,7 @@ #ifndef TIMER_H #define TIMER_H -#include "stdio.h" -#include "stdlib.h" -#include "time.h" +#include class Timer { public: diff --git a/tools/phonon/tricubic/CMakeLists.txt b/tools/phonon/tricubic/CMakeLists.txt new file mode 100644 index 0000000000..190c6ceb80 --- /dev/null +++ b/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}) diff --git a/tools/phonon/tricubic/LICENSE b/tools/phonon/tricubic/LICENSE new file mode 100644 index 0000000000..d60c31a97a --- /dev/null +++ b/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. + + + Copyright (C) + + 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. + + , 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. diff --git a/tools/phonon/tricubic/README.md b/tools/phonon/tricubic/README.md new file mode 100644 index 0000000000..0e8058339e --- /dev/null +++ b/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. + diff --git a/tools/phonon/tricubic/tricubic.cpp b/tools/phonon/tricubic/tricubic.cpp new file mode 100644 index 0000000000..974ed6aa52 --- /dev/null +++ b/tools/phonon/tricubic/tricubic.cpp @@ -0,0 +1,185 @@ + +#include "tricubic.h" + +#include + +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