Added unit test.

materialsproject · Jul 2, 2018 · 9cbebb8 · 9cbebb8
1 parent 82380ca
commit 9cbebb8
Show file tree

Hide file tree

Showing 2 changed files with 28 additions and 7 deletions.
diff --git a/pymatgen/analysis/local_env.py b/pymatgen/analysis/local_env.py
@@ -2376,15 +2376,18 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
                         tmp = max(
                             -1.0, min(np.inner(zaxis, rijnorm[k]), 1.0))
                         thetak = acos(tmp)
-                        xaxistmp = gramschmidt(rijnorm[k], zaxis)
-                        if np.linalg.norm(xaxistmp) < very_small:
+                        xaxis = gramschmidt(rijnorm[k], zaxis)
+                        if np.linalg.norm(xaxis) < very_small:
                             flag_xaxis = True
                         else:
-                            xaxis = xaxistmp / np.linalg.norm(xaxistmp)
+                            xaxis = xaxis / np.linalg.norm(xaxis)
                             flag_xaxis = False
                         if self._comp_azi:
-                            yaxistmp = np.cross(zaxis, xaxis)
-                            yaxis = yaxistmp / np.linalg.norm(yaxistmp)
+                            flag_yaxis = True
+                            yaxis = np.cross(zaxis, xaxis)
+                            if np.linalg.norm(yaxis) > very_small:
+                                yaxis = yaxis / np.linalg.norm(yaxis)
+                                flag_yaxis = False
 
                         # Contributions of j-i-k angles, where i represents the
                         # central atom and j and k two of the neighbors.
@@ -2589,7 +2592,7 @@ def get_order_parameters(self, structure, n, indices_neighs=None, \
                                                     qsptheta[i][j][kc] += exp(-0.5 * tmp * tmp) * \
                                                             exp(-0.5 * tmp2 * tmp2)
                                                     norms[i][j][kc] += 1.0
-                                        elif t == "sq_face_cap_trig_pris":
+                                        elif t == "sq_face_cap_trig_pris" and not flag_yaxis:
                                             if thetak < self._params[i]['TA3']:
                                                 if thetam < self._params[i]['TA3']:
                                                     tmp = cos(self._params[i]['fac_AA1'] * \

diff --git a/pymatgen/analysis/tests/test_local_env.py b/pymatgen/analysis/tests/test_local_env.py
@@ -696,6 +696,18 @@ def setUp(self):
             [0.0, 0.0, -1.0], [-1.0, 0.0, 0.0]],
             validate_proximity=False, to_unit_cell=False,
             coords_are_cartesian=True, site_properties=None)
+        self.sq_face_capped_trig_pris = Structure(
+            Lattice.from_lengths_and_angles(
+            [30, 30, 30], [90, 90, 90]),
+            ["H", "H", "H", "H", "H", "H", "H", "H"],
+            [[0, 0, 0], [-0.6546536707079771, -0.37796447300922725, 0.6546536707079771],
+            [0.6546536707079771, -0.37796447300922725, 0.6546536707079771],
+            [0.0, 0.7559289460184545, 0.6546536707079771],
+            [-0.6546536707079771, -0.37796447300922725, -0.6546536707079771],
+            [0.6546536707079771, -0.37796447300922725, -0.6546536707079771],
+            [0.0, 0.7559289460184545, -0.6546536707079771], [0.0, -1.0, 0.0]],
+            validate_proximity=False, to_unit_cell=False,
+            coords_are_cartesian=True, site_properties=None)
 
     def test_init(self):
         self.assertIsNotNone(
@@ -714,7 +726,7 @@ def test_get_order_parameters(self):
             "tri_plan", "sq_plan", "pent_plan", "sq_pyr", "tri_pyr", \
             "pent_pyr", "hex_pyr", "pent_bipyr", "hex_bipyr", "T", "cuboct", \
             "see_saw_rect", "hex_plan_max", "tet_max", "oct_max", "tri_plan_max", "sq_plan_max", \
-            "pent_plan_max", "cuboct_max", "tet_max"]
+            "pent_plan_max", "cuboct_max", "tet_max", "sq_face_cap_trig_pris"]
         op_params = [None for i in range(len(op_types))]
         op_params[1] = {'TA': 1, 'IGW_TA': 1./0.0667}
         op_params[2] = {'TA': 45./180, 'IGW_TA': 1./0.0667}
@@ -893,6 +905,12 @@ def test_get_order_parameters(self):
             self.hexagonal_planar, 0, indices_neighs=[1,2,3,4,5,6])
         self.assertAlmostEqual(int(op_vals[26] * 1000 + 0.5), 1000)
 
+        # Square face capped trigonal prism.
+        op_vals = ops_101.get_order_parameters(
+            self.sq_face_capped_trig_pris, 0,
+            indices_neighs=[i for i in range(1, 8)])
+        self.assertAlmostEqual(int(op_vals[34] * 1000 + 0.5), 1000)
+
         # Test providing explicit neighbor lists.
         op_vals = ops_101.get_order_parameters(self.bcc, 0, indices_neighs=[1])
         self.assertIsNotNone(op_vals[0])