From 01fcf30ec0a53165480c4f7117a2d90b55864a0a Mon Sep 17 00:00:00 2001 From: bpuchala Date: Thu, 17 Aug 2023 09:21:48 -0400 Subject: [PATCH] fix abs and pow usage --- include/casm/mapping/impl/LatticeMap.hh | 6 +-- python/tests/test_map_structures.py | 50 ++++++++++++++++++++++++- src/casm/mapping/atom_cost.cc | 12 +++--- src/casm/mapping/impl/LatticeMap.cc | 2 +- src/casm/mapping/impl/StrucMapping.cc | 5 ++- src/casm/mapping/lattice_cost.cc | 2 +- 6 files changed, 64 insertions(+), 13 deletions(-) diff --git a/include/casm/mapping/impl/LatticeMap.hh b/include/casm/mapping/impl/LatticeMap.hh index f6b01ba..30b2178 100644 --- a/include/casm/mapping/impl/LatticeMap.hh +++ b/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 diff --git a/python/tests/test_map_structures.py b/python/tests/test_map_structures.py index 4f65851..cf5371e 100644 --- a/python/tests/test_map_structures.py +++ b/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) diff --git a/src/casm/mapping/atom_cost.cc b/src/casm/mapping/atom_cost.cc index 2dfc7f3..29a5e94 100644 --- a/src/casm/mapping/atom_cost.cc +++ b/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 diff --git a/src/casm/mapping/impl/LatticeMap.cc b/src/casm/mapping/impl/LatticeMap.cc index a5defc8..804082d 100644 --- a/src/casm/mapping/impl/LatticeMap.cc +++ b/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 = diff --git a/src/casm/mapping/impl/StrucMapping.cc b/src/casm/mapping/impl/StrucMapping.cc index e01e0fb..faca3b0 100644 --- a/src/casm/mapping/impl/StrucMapping.cc +++ b/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); } diff --git a/src/casm/mapping/lattice_cost.cc b/src/casm/mapping/lattice_cost.cc index f7d6f71..6bdcd51 100644 --- a/src/casm/mapping/lattice_cost.cc +++ b/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);