Add test for get_fc3 with lapacke pinv solver

phono3py/phonon3/fc3.py

@@ -65,14 +65,20 @@ def get_fc3(
     disp_dataset,
     symmetry: Symmetry,
     is_compact_fc=False,
+    pinv_solver: str = "numpy",
     verbose=False,
 ):
     """Calculate fc3.
 
-    Even when 'cutoff_distance' in dataset, all displacements are in the dataset,
-    but force-sets out of cutoff-pair-distance are zero. fc3 is solved in exactly
-    the same way. Then post-clean-up
-    is performed.
+    Even when 'cutoff_distance' in dataset, all displacements are in the
+    dataset, but force-sets out of cutoff-pair-distance are zero. fc3 is solved
+    in exactly the same way. Then post-clean-up is performed.
+
+    Returns
+    -------
+    tuple :
+        (fc2, fc3) fc2 and fc3 can be compact or full array formats depending on
+        `is_compact_fc`. See Phono3py.produce_fc3.
 
     """
     # fc2 has to be full matrix to compute delta-fc2
@@ -85,6 +91,7 @@ def get_fc3(
         fc2,
         symmetry,
         is_compact_fc=(is_compact_fc and "cutoff_distance" not in disp_dataset),
+        pinv_solver=pinv_solver,
         verbose=verbose,
     )
     if verbose:
@@ -669,6 +676,7 @@ def _get_fc3_least_atoms(
     fc2,
     symmetry: Symmetry,
     is_compact_fc: bool = False,
+    pinv_solver="numpy",
     verbose: bool = True,
 ):
     symprec = symmetry.tolerance
@@ -735,6 +743,7 @@ def _get_fc3_least_atoms(
             displacements_first,
             np.array(delta_fc2s, dtype="double", order="C"),
             symprec,
+            pinv_solver=pinv_solver,
             verbose=verbose,
         )
         if is_compact_fc:

test/phonon3/test_fc3.py

@@ -1,8 +1,14 @@
 """Tests for fc3."""
 import numpy as np
+import pytest
 
 from phono3py import Phono3py
-from phono3py.phonon3.fc3 import cutoff_fc3_by_zero
+from phono3py.phonon3.fc3 import (
+    cutoff_fc3_by_zero,
+    get_fc3,
+    set_permutation_symmetry_fc3,
+    set_translational_invariance_fc3,
+)
 
 
 def test_cutoff_fc3(nacl_pbe_cutoff_fc3: Phono3py, nacl_pbe: Phono3py):
@@ -90,6 +96,42 @@ def test_fc3(si_pbesol_111: Phono3py):
     np.testing.assert_allclose(ph.fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0)
 
 
+@pytest.mark.parametrize("pinv_solver", ["numpy", "lapacke"])
+def test_fc3_lapacke_solver(si_pbesol_111: Phono3py, pinv_solver: str):
+    """Test fc3 with Si PBEsol 1x1x1 using lapacke solver."""
+    ph = si_pbesol_111
+    _, fc3 = get_fc3(
+        ph.supercell,
+        ph.primitive,
+        ph.dataset,
+        ph.symmetry,
+        pinv_solver=pinv_solver,
+        verbose=True,
+    )
+    set_translational_invariance_fc3(fc3)
+    set_permutation_symmetry_fc3(fc3)
+
+    fc3_ref = [
+        [
+            [1.07250822e-01, 1.86302073e-17, -4.26452855e-18],
+            [8.96414569e-03, -1.43046911e-01, -1.38498937e-01],
+            [-8.96414569e-03, -1.38498937e-01, -1.43046911e-01],
+        ],
+        [
+            [-8.96414569e-03, -1.43046911e-01, -1.38498937e-01],
+            [-3.39457157e-02, -4.63315728e-17, -4.17779237e-17],
+            [-3.31746167e-01, -2.60025724e-02, -2.60025724e-02],
+        ],
+        [
+            [8.96414569e-03, -1.38498937e-01, -1.43046911e-01],
+            [-3.31746167e-01, 2.60025724e-02, 2.60025724e-02],
+            [-3.39457157e-02, 3.69351540e-17, 5.94504191e-18],
+        ],
+    ]
+
+    np.testing.assert_allclose(fc3[0, 1, 7], fc3_ref, atol=1e-8, rtol=0)
+
+
 # @pytest.mark.skipif(not FC_CALCULATOR_ALM_AVAILABLE, reason="not found ALM package")
 def test_fc3_alm(si_pbesol_111_alm: Phono3py):
     """Test fc3 with Si PBEsol 1x1x1 calcualted using ALM."""