fix abs and pow usage

include/casm/mapping/impl/LatticeMap.hh

@@ -31,10 +31,10 @@ class StrainCostCalculator {
       Eigen::Matrix3d const &_deformation_gradient, double _vol_factor);
 
   // \brief Volumetric factor :
-  // pow(abs(_deformation_gradient.determinant()),1./3.), used to normalize the
-  // strain cost to make it volume-independent
+  // std::pow(std::abs(_deformation_gradient.determinant()),1./3.), used to
+  // normalize the strain cost to make it volume-independent
   static double vol_factor(Eigen::Matrix3d const &_deformation_gradient) {
-    return pow(std::abs(_deformation_gradient.determinant()), 1. / 3.);
+    return std::pow(std::abs(_deformation_gradient.determinant()), 1. / 3.);
   }
 
   //\brief Anisotropic strain cost; utilizes stored gram matrix to compute

python/tests/test_map_structures.py

@@ -8,6 +8,49 @@
 import libcasm.xtal.structures as xtal_structures
 
 
+def make_geometric_atom_cost(
+    supercell_lattice_column_vector_matrix,
+    displacement,
+):
+    S = supercell_lattice_column_vector_matrix
+    N_site = displacement.shape[1]
+    volume_per_site = np.abs(np.linalg.det(S)) / N_site
+    displacement_squaredNorm = np.sum(displacement**2)
+    geometric_atom_cost = (
+        math.pow(3 * volume_per_site / (4.0 * np.pi), -2.0 / 3.0)
+        * displacement_squaredNorm
+        / N_site
+    )
+    # print("volume_per_site:", volume_per_site)
+    # print("displacement_squaredNorm:", displacement_squaredNorm)
+    return geometric_atom_cost
+
+
+def make_isotropic_atom_cost(
+    prim_lattice_column_vector_matrix,
+    lattice_mapping,
+    displacement,
+):
+    L1 = prim_lattice_column_vector_matrix
+    T = lattice_mapping.transformation_matrix_to_super()
+    N = lattice_mapping.reorientation()
+    U = lattice_mapping.right_stretch()
+
+    S1 = L1 @ T @ N
+    L2 = U @ S1
+    d = displacement
+    d_reverse = -U @ d
+
+    geometric_atom_cost_forward = make_geometric_atom_cost(S1, d)
+    geometric_atom_cost_reverse = make_geometric_atom_cost(L2, d_reverse)
+    isotropic_atom_cost = (
+        geometric_atom_cost_forward + geometric_atom_cost_reverse
+    ) / 2.0
+    # print("geometric_atom_cost_forward:", geometric_atom_cost_forward)
+    # print("geometric_atom_cost_reverse:", geometric_atom_cost_reverse)
+    return isotropic_atom_cost
+
+
 def check_mapping(prim, structure, structure_mapping):
     # print("structure:")
     # print("lattice_column_vector_matrix:\n",
@@ -76,6 +119,11 @@ def check_mapping(prim, structure, structure_mapping):
             x2 = r2[:, perm[i]] + trans
             d = xtal.min_periodic_displacement(structure.lattice(), x1, x2)
             assert math.isclose(np.linalg.norm(d), 0.0, abs_tol=1e-10)
+
+    # atom mapping cost:
+    # isotropic_atom_cost = make_isotropic_atom_cost(L1, lmap, disp)
+    # print("isotropic_atom_cost:", isotropic_atom_cost)
+
     # assert True == False
 
 
@@ -339,6 +387,6 @@ def test_bcc_hcp_mapping():
 
     assert len(structure_mappings)
     for i, smap in enumerate(structure_mappings):
+        check_mapping(prim, hcp_structure, smap)
         assert math.isclose(smap.lattice_cost(), 0.007297079413597657)
         assert math.isclose(smap.atom_cost(), 0.06274848406141671)
-        check_mapping(prim, hcp_structure, smap)

src/casm/mapping/atom_cost.cc

@@ -29,8 +29,8 @@ double make_geometric_atom_cost(
     Eigen::MatrixXd const &displacement) {
   double N_site = displacement.cols();
   double volume_per_site =
-      abs(supercell_lattice_column_vector_matrix.determinant()) / N_site;
-  return pow(3. * volume_per_site / (4. * M_PI), -2. / 3.) *
+      std::abs(supercell_lattice_column_vector_matrix.determinant()) / N_site;
+  return std::pow(3. * volume_per_site / (4. * M_PI), -2. / 3.) *
          displacement.squaredNorm() / N_site;
 }
 
@@ -89,9 +89,11 @@ double make_isotropic_atom_cost(
   Eigen::MatrixXd const &d = displacement;
   Eigen::MatrixXd d_reverse = -U * d;
 
-  return (make_geometric_atom_cost(S1, d) +
-          make_geometric_atom_cost(L2, d_reverse)) /
-         2.;
+  double isotropic_atom_cost = (make_geometric_atom_cost(S1, d) +
+                                make_geometric_atom_cost(L2, d_reverse)) /
+                               2.;
+  std::cout << "cxx isotropic_atom_cost: " << isotropic_atom_cost << std::endl;
+  return isotropic_atom_cost;
 }
 
 /// \brief Return the symmetry-preserving component of displacement

src/casm/mapping/impl/LatticeMap.cc

@@ -686,7 +686,7 @@ bool LatticeMap::_check_canonical() const {
 double isotropic_strain_cost(Eigen::Matrix3d const &deformation_gradient) {
   // written using convention B = V - I of the mapping paper:
   Eigen::Matrix3d const &F_reverse = deformation_gradient;
-  double vol_factor = pow(std::abs(F_reverse.determinant()), 1. / 3.);
+  double vol_factor = std::pow(std::abs(F_reverse.determinant()), 1. / 3.);
 
   Eigen::Matrix3d I = Eigen::Matrix3d::Identity(3, 3);
   Eigen::Matrix3d U_reverse_normalized =

src/casm/mapping/impl/StrucMapping.cc

@@ -35,7 +35,7 @@ double atomic_cost_child(const MappingNode &mapped_result, Index Nsites) {
   double atomic_vol =
       mapped_result.lattice_node.parent.superlattice().volume() /
       double(Nsites) / mapped_result.lattice_node.stretch.determinant();
-  return pow(3. * abs(atomic_vol) / (4. * M_PI), -2. / 3.) *
+  return std::pow(3. * std::abs(atomic_vol) / (4. * M_PI), -2. / 3.) *
          (mapped_result.lattice_node.stretch.inverse() *
           mapped_result.atom_displacement)
              .squaredNorm() /
@@ -49,7 +49,8 @@ double atomic_cost_parent(const MappingNode &mapped_result, Index Nsites) {
   double atomic_vol =
       mapped_result.lattice_node.parent.superlattice().volume() /
       double(Nsites);
-  return pow(3. * abs(atomic_vol) / (4. * M_PI), -2. / 3.) *
+
+  return std::pow(3. * std::abs(atomic_vol) / (4. * M_PI), -2. / 3.) *
          (mapped_result.atom_displacement).squaredNorm() / double(Nsites);
 }
 

src/casm/mapping/lattice_cost.cc

@@ -51,7 +51,7 @@ namespace mapping {
 ///
 double isotropic_strain_cost(Eigen::Matrix3d const &deformation_gradient) {
   Eigen::Matrix3d const &F = deformation_gradient;
-  double vol_factor = pow(std::abs(F.determinant()), 1. / 3.);
+  double vol_factor = std::pow(std::abs(F.determinant()), 1. / 3.);
 
   Eigen::Matrix3d I = Eigen::Matrix3d::Identity(3, 3);
   Eigen::Matrix3d U_normalized = polar_decomposition(F / vol_factor);