Fix ddm with Wang's NAC method

c/_phonopy.cpp

@@ -247,15 +247,17 @@ void py_get_recip_dipole_dipole_q0(nb::ndarray<> py_dd_q0,
 
 void py_get_derivative_dynmat(
     nb::ndarray<> py_derivative_dynmat, nb::ndarray<> py_force_constants,
-    nb::ndarray<> py_q_vector, nb::ndarray<> py_lattice, nb::ndarray<> py_svecs,
+    nb::ndarray<> py_q_vector, nb::ndarray<> py_lattice,
-    nb::ndarray<> py_multi, nb::ndarray<> py_masses, nb::ndarray<> py_s2p_map,
+    nb::ndarray<> py_reclat, nb::ndarray<> py_svecs, nb::ndarray<> py_multi,
-    nb::ndarray<> py_p2s_map, double nac_factor, nb::ndarray<> py_born,
+    nb::ndarray<> py_masses, nb::ndarray<> py_s2p_map, nb::ndarray<> py_p2s_map,
-    nb::ndarray<> py_dielectric, nb::ndarray<> py_q_direction, long is_nac,
+    double nac_factor, nb::ndarray<> py_born, nb::ndarray<> py_dielectric,
-    long is_nac_q_zero, long use_openmp) {
+    nb::ndarray<> py_q_direction, long is_nac, long is_nac_q_zero,
+    long use_openmp) {
     double(*ddm)[2];
     double *fc;
     double *q_vector;
-    double *lat;
+    double *lattice;
+    double *reclat;
     double(*svecs)[3];
     double *masses;
     long(*multi)[2];
@@ -271,7 +273,8 @@ void py_get_derivative_dynmat(
     ddm = (double(*)[2])py_derivative_dynmat.data();
     fc = (double *)py_force_constants.data();
     q_vector = (double *)py_q_vector.data();
-    lat = (double *)py_lattice.data();
+    lattice = (double *)py_lattice.data();
+    reclat = (double *)py_reclat.data();
     svecs = (double(*)[3])py_svecs.data();
     masses = (double *)py_masses.data();
     multi = (long(*)[2])py_multi.data();
@@ -288,9 +291,10 @@ void py_get_derivative_dynmat(
         q_dir = (double *)py_q_direction.data();
     }
 
-    phpy_get_derivative_dynmat_at_q(
+    phpy_get_derivative_dynmat_at_q(ddm, num_patom, num_satom, fc, q_vector,
-        ddm, num_patom, num_satom, fc, q_vector, lat, svecs, multi, masses,
+                                    lattice, reclat, svecs, multi, masses,
-        s2p_map, p2s_map, nac_factor, born, epsilon, q_dir, is_nac, use_openmp);
+                                    s2p_map, p2s_map, nac_factor, born, epsilon,
+                                    q_dir, is_nac, use_openmp);
 }
 
 void py_get_thermal_properties(nb::ndarray<> py_thermal_props,

c/derivative_dynmat.c

@@ -35,6 +35,7 @@
 #include "derivative_dynmat.h"
 
 #include <math.h>
+#include <stdio.h>
 #include <stdlib.h>
 #define PI 3.14159265358979323846
 
@@ -44,8 +45,7 @@ static void get_derivative_dynmat_at_q(
     const long num_patom, const long num_satom, const double *fc,
     const double *q, const double *lattice, /* column vector */
     const double (*svecs)[3], const long (*multi)[2], const double *mass,
-    const long *s2p_map, const long *p2s_map, const double nac_factor,
+    const long *s2p_map, const long *p2s_map);
-    const double *born, const double *dielectric);
 static void get_derivative_nac(double *ddnac, double *dnac,
                                const long num_patom, const double *lattice,
                                const double *mass, const double *q,
@@ -63,6 +63,7 @@ void ddm_get_derivative_dynmat_at_q(
     double (*derivative_dynmat)[2], const long num_patom, const long num_satom,
     const double *fc, const double *q,
     const double *lattice, /* column vector */
+    const double *reclat,  /* column vector */
     const double (*svecs)[3], const long (*multi)[2], const double *mass,
     const long *s2p_map, const long *p2s_map, const double nac_factor,
     const double *born, const double *dielectric, const double *q_direction,
@@ -74,8 +75,8 @@ void ddm_get_derivative_dynmat_at_q(
     if (is_nac) {
         ddnac = (double *)malloc(sizeof(double) * num_patom * num_patom * 27);
         dnac = (double *)malloc(sizeof(double) * num_patom * num_patom * 9);
-        factor = nac_factor * num_patom / num_satom;
+        factor = (nac_factor * num_patom) / num_satom;
-        get_derivative_nac(ddnac, dnac, num_patom, lattice, mass, q, born,
+        get_derivative_nac(ddnac, dnac, num_patom, reclat, mass, q, born,
                            dielectric, q_direction, factor);
     }
 
@@ -89,15 +90,15 @@ void ddm_get_derivative_dynmat_at_q(
             get_derivative_dynmat_at_q(derivative_dynmat, i, j, ddnac, dnac,
                                        is_nac, num_patom, num_satom, fc, q,
                                        lattice, svecs, multi, mass, s2p_map,
-                                       p2s_map, nac_factor, born, dielectric);
+                                       p2s_map);
         }
     } else {
         for (i = 0; i < num_patom; i++) {
             for (j = 0; j < num_patom; j++) {
-                get_derivative_dynmat_at_q(
+                get_derivative_dynmat_at_q(derivative_dynmat, i, j, ddnac, dnac,
-                    derivative_dynmat, i, j, ddnac, dnac, is_nac, num_patom,
+                                           is_nac, num_patom, num_satom, fc, q,
-                    num_satom, fc, q, lattice, svecs, multi, mass, s2p_map,
+                                           lattice, svecs, multi, mass, s2p_map,
-                    p2s_map, nac_factor, born, dielectric);
+                                           p2s_map);
             }
         }
     }
@@ -126,14 +127,15 @@ void ddm_get_derivative_dynmat_at_q(
     }
 }
 
-void get_derivative_dynmat_at_q(
+void get_derivative_dynmat_at_q(double (*derivative_dynmat)[2], const long i,
-    double (*derivative_dynmat)[2], const long i, const long j,
+                                const long j, const double *ddnac,
-    const double *ddnac, const double *dnac, const long is_nac,
+                                const double *dnac, const long is_nac,
-    const long num_patom, const long num_satom, const double *fc,
+                                const long num_patom, const long num_satom,
-    const double *q, const double *lattice, /* column vector */
+                                const double *fc, const double *q,
-    const double (*svecs)[3], const long (*multi)[2], const double *mass,
+                                const double *lattice, /* column vector */
-    const long *s2p_map, const long *p2s_map, const double nac_factor,
+                                const double (*svecs)[3],
-    const double *born, const double *dielectric) {
+                                const long (*multi)[2], const double *mass,
+                                const long *s2p_map, const long *p2s_map) {
     long k, l, m, n, adrs, m_pair, i_pair, svecs_adrs;
     double coef[3], real_coef[3], imag_coef[3];
     double c, s, phase, mass_sqrt, fc_elem, real_phase, imag_phase;
@@ -237,35 +239,21 @@ void get_derivative_dynmat_at_q(
 /* D_nac = a * AB/C */
 /* dD_nac = a * D_nac * (A'/A + B'/B - C'/C) */
 static void get_derivative_nac(double *ddnac, double *dnac,
-                               const long num_patom, const double *lattice,
+                               const long num_patom, const double *reclat,
                                const double *mass, const double *q,
                                const double *born, const double *dielectric,
                                const double *q_direction, const double factor) {
     long i, j, k, l, m;
-    double a, b, c, da, db, dc, volume, mass_sqrt;
+    double a, b, c, da, db, dc, mass_sqrt;
-    double q_cart[3], rec_lat[9];
+    double q_cart[3];
-
-    volume = lattice[0] * (lattice[4] * lattice[8] - lattice[5] * lattice[7]) +
-             lattice[1] * (lattice[5] * lattice[6] - lattice[3] * lattice[8]) +
-             lattice[2] * (lattice[3] * lattice[7] - lattice[4] * lattice[6]);
-
-    rec_lat[0] = (lattice[4] * lattice[8] - lattice[5] * lattice[7]) / volume;
-    rec_lat[1] = (lattice[5] * lattice[6] - lattice[3] * lattice[8]) / volume;
-    rec_lat[2] = (lattice[3] * lattice[7] - lattice[4] * lattice[6]) / volume;
-    rec_lat[3] = (lattice[7] * lattice[2] - lattice[8] * lattice[1]) / volume;
-    rec_lat[4] = (lattice[8] * lattice[0] - lattice[6] * lattice[2]) / volume;
-    rec_lat[5] = (lattice[6] * lattice[1] - lattice[7] * lattice[0]) / volume;
-    rec_lat[6] = (lattice[1] * lattice[5] - lattice[2] * lattice[4]) / volume;
-    rec_lat[7] = (lattice[2] * lattice[3] - lattice[0] * lattice[5]) / volume;
-    rec_lat[8] = (lattice[0] * lattice[4] - lattice[1] * lattice[3]) / volume;
 
     for (i = 0; i < 3; i++) {
         q_cart[i] = 0;
         for (j = 0; j < 3; j++) {
             if (q_direction) {
-                q_cart[i] += rec_lat[i * 3 + j] * q_direction[j];
+                q_cart[i] += reclat[i * 3 + j] * q_direction[j];
             } else {
-                q_cart[i] += rec_lat[i * 3 + j] * q[j];
+                q_cart[i] += reclat[i * 3 + j] * q[j];
             }
         }
     }

c/derivative_dynmat.h

@@ -39,6 +39,7 @@ void ddm_get_derivative_dynmat_at_q(
     double (*derivative_dynmat)[2], const long num_patom, const long num_satom,
     const double *fc, const double *q,
     const double *lattice, /* column vector */
+    const double *reclat,  /* column vector */
     const double (*svecs)[3], const long (*multi)[2], const double *mass,
     const long *s2p_map, const long *p2s_map, const double nac_factor,
     const double *born, const double *dielectric, const double *q_direction,

c/phonopy.c

@@ -134,14 +134,15 @@ void phpy_get_derivative_dynmat_at_q(
     double (*derivative_dynmat)[2], const long num_patom, const long num_satom,
     const double *fc, const double *q,
     const double *lattice, /* column vector */
+    const double *reclat,  /* column vector */
     const double (*svecs)[3], const long (*multi)[2], const double *mass,
     const long *s2p_map, const long *p2s_map, const double nac_factor,
     const double *born, const double *dielectric, const double *q_direction,
     const long is_nac, const long use_openmp) {
     ddm_get_derivative_dynmat_at_q(derivative_dynmat, num_patom, num_satom, fc,
-                                   q, lattice, svecs, multi, mass, s2p_map,
+                                   q, lattice, reclat, svecs, multi, mass,
-                                   p2s_map, nac_factor, born, dielectric,
+                                   s2p_map, p2s_map, nac_factor, born,
-                                   q_direction, is_nac, use_openmp);
+                                   dielectric, q_direction, is_nac, use_openmp);
 }
 
 void phpy_get_relative_grid_address(long relative_grid_address[24][4][3],

c/phonopy.h

@@ -82,6 +82,7 @@ void phpy_get_derivative_dynmat_at_q(
     double (*derivative_dynmat)[2], const long num_patom, const long num_satom,
     const double *fc, const double *q,
     const double *lattice, /* column vector */
+    const double *reclat,  /* column vector */
     const double (*svecs)[3], const long (*multi)[2], const double *mass,
     const long *s2p_map, const long *p2s_map, const double nac_factor,
     const double *born, const double *dielectric, const double *q_direction,

phonopy/harmonic/derivative_dynmat.py

@@ -55,6 +55,8 @@ class DerivativeOfDynamicalMatrix:
 
     """
 
+    Q_DIRECTION_TOLERANCE = 1e-5
+
     def __init__(self, dynamical_matrix: Union[DynamicalMatrix, DynamicalMatrixNAC]):
         """Init method.
 
@@ -83,7 +85,6 @@ class DerivativeOfDynamicalMatrix:
             self._multi = multi
         else:
             self._svecs, self._multi = sparse_to_dense_svecs(svecs, multi)
-        # self._svecs, self._multi = self._pcell.get_smallest_vectors()
 
         self._ddm = None
 
@@ -134,30 +135,34 @@ class DerivativeOfDynamicalMatrix:
         import phonopy._phonopy as phonoc
 
         num_patom = len(self._p2s_map)
-
         fc = self._force_constants
         ddm = np.zeros(
             (3, num_patom * 3, num_patom * 3),
             dtype=("c%d" % (np.dtype("double").itemsize * 2)),
         )
+        reclat = np.array(np.linalg.inv(self._pcell.cell), dtype="double", order="C")
+
+        is_nac = False
+        is_nac_q_zero = True
+        born = np.zeros(9)  # dummy value
+        dielectric = np.zeros(9)  # dummy value
+        nac_factor = 0  # dummy value
+        q_dir = np.zeros(3)  # dummy value
+
         if isinstance(self._dynmat, DynamicalMatrixNAC):
             is_nac = True
             born = self._dynmat.born
             dielectric = self._dynmat.dielectric_constant
             nac_factor = self._dynmat.nac_factor
             if q_direction is None:
-                q_dir = np.zeros(3)
+                q_norm = np.linalg.norm(reclat @ q)
-                is_nac_q_zero = True
+                if q_norm < self.Q_DIRECTION_TOLERANCE:
-            else:
-                q_dir = np.array(q_direction, dtype="double", order="C")
-                is_nac_q_zero = False
-        else:
-            born = np.zeros(0)  # dummy value
-            dielectric = np.zeros(0)  # dummy value
-            nac_factor = 0  # dummy value
-            q_dir = np.zeros(3)  # dummy value
                     is_nac = False
+                else:
                     is_nac_q_zero = True
+            else:
+                q_dir = np.array(q_direction, dtype="double")
+                is_nac_q_zero = False
 
         if fc.shape[0] == fc.shape[1]:  # full fc
             phonoc.derivative_dynmat(
@@ -165,6 +170,7 @@ class DerivativeOfDynamicalMatrix:
                 fc,
                 np.array(q, dtype="double"),
                 np.array(self._pcell.cell.T, dtype="double", order="C"),
+                reclat,
                 self._svecs,
                 self._multi,
                 self._pcell.masses,
@@ -184,6 +190,7 @@ class DerivativeOfDynamicalMatrix:
                 fc,
                 np.array(q, dtype="double"),
                 np.array(self._pcell.cell.T, dtype="double", order="C"),
+                np.array(np.linalg.inv(self._pcell.cell), dtype="double", order="C"),
                 self._svecs,
                 self._multi,
                 self._pcell.masses,
@@ -236,7 +243,7 @@ class DerivativeOfDynamicalMatrix:
                     continue
 
                 multi = multiplicity[k, i]
-                vecs_multi = vecs[k, i, :multi]
+                vecs_multi = vecs[multi[1] : multi[1] + multi[0]]
                 phase_multi = np.exp(
                     [np.vdot(vec, q) * 2j * np.pi for vec in vecs_multi]
                 )
@@ -263,7 +270,7 @@ class DerivativeOfDynamicalMatrix:
                     ):
                         ddm_elem += d_nac[ll, i, j] * phase_multi.sum()
 
-                    ddm_local[ll] += ddm_elem / mass / multi
+                    ddm_local[ll] += ddm_elem / mass / multi[0]
 
             ddm[:, (i * 3) : (i * 3 + 3), (j * 3) : (j * 3 + 3)] = ddm_local
 

test/conftest.py

@@ -43,6 +43,26 @@ def ph_nacl() -> Phonopy:
     )
 
 
+@pytest.fixture(scope="session")
+def ph_nacl_wang() -> Phonopy:
+    """Return Phonopy class instance of NaCl 2x2x2."""
+    yaml_filename = cwd / "phonopy_disp_NaCl.yaml"
+    force_sets_filename = cwd / "FORCE_SETS_NaCl"
+    born_filename = cwd / "BORN_NaCl"
+    ph = phonopy.load(
+        yaml_filename,
+        force_sets_filename=force_sets_filename,
+        born_filename=born_filename,
+        is_compact_fc=False,
+        log_level=1,
+        produce_fc=True,
+    )
+    nac_params = ph.nac_params
+    nac_params["method"] = "wang"
+    ph.nac_params = nac_params
+    return ph
+
+
 @pytest.fixture(scope="session")
 def ph_nacl_nofcsym() -> Phonopy:
     """Return Phonopy class instance of NaCl 2x2x2 without symmetrizing fc2."""

test/harmonic/test_derivative_dynmat.py

@@ -6,6 +6,7 @@ from phonopy import Phonopy
 from phonopy.harmonic.derivative_dynmat import DerivativeOfDynamicalMatrix
 
 ddm_ph_nacl = [
+    [
         -0.2568036,
         0.6947168,
         0.0337825,
@@ -34,9 +35,41 @@ ddm_ph_nacl = [
         1.3859454,
         -1.5152406,
         1.3859454,
+    ],
+    [
+        -0.4582279,
+        1.1733186,
+        0.0783508,
+        0.0543731,
+        0.0783508,
+        0.0543731,
+        1.1733186,
+        0.0724132,
+        0.0543731,
+        0.2215595,
+        0.0543731,
+        0.2215595,
+        0.5943896,
+        -0.2136538,
+        0.5943896,
+        -0.2136538,
+        -0.2136538,
+        0.4540728,
+        -0.2136538,
+        0.4540728,
+        0.5943896,
+        -0.2136538,
+        0.5943896,
+        -0.2136538,
+        -0.2136538,
+        0.4540728,
+        -0.2136538,
+        0.4540728,
+    ],
 ]
 
 ddm_ph_nacl_nonac = [
+    [
         0.8832905,
         -0.2233653,
         0.0337825,
@@ -65,6 +98,37 @@ ddm_ph_nacl_nonac = [
         0.2668973,
         -0.1255825,
         0.2668973,
+    ],
+    [
+        0.2826501,
+        0.5767127,
+        0.0783508,
+        0.0543731,
+        0.0783508,
+        0.0543731,
+        0.5767127,
+        0.5528413,
+        0.0543731,
+        0.2215595,
+        0.0543731,
+        0.2215595,
+        0.5943896,
+        -0.2136538,
+        0.5943896,
+        -0.2136538,
+        -0.2136538,
+        0.4540728,
+        -0.2136538,
+        0.4540728,
+        0.5943896,
+        -0.2136538,
+        0.5943896,
+        -0.2136538,
+        -0.2136538,
+        0.4540728,
+        -0.2136538,
+        0.4540728,
+    ],
 ]
 
 
@@ -73,6 +137,11 @@ def test_ddm_nac(ph_nacl: Phonopy):
     _assert(ph_nacl, ddm_ph_nacl)
 
 
+def test_ddm_nac_wang(ph_nacl_wang: Phonopy):
+    """Test by NaCl with Wang's NAC."""
+    _assert(ph_nacl_wang, ddm_ph_nacl, show=True)
+
+
 def test_ddm_nac_compact(ph_nacl_compact_fcsym: Phonopy):
     """Test by NaCl with compact fc."""
     _assert(ph_nacl_compact_fcsym, ddm_ph_nacl)
@@ -91,13 +160,14 @@ def test_ddm_nonac_compact(ph_nacl_nonac_compact_fc: Phonopy):
 def _assert(ph: Phonopy, ref_vals, show=False):
     dynmat = ph.dynamical_matrix
     ddynmat = DerivativeOfDynamicalMatrix(dynmat)
-    ddynmat.run([0, 0.1, 0.1])
+    for i, q in enumerate(([0, 0.1, 0.1], [0, 0.25, 0.25])):
+        ddynmat.run(q)
         ddm = ddynmat.d_dynamical_matrix
         condition = np.abs(ddm) > 1e-8
         vals = np.extract(condition, ddm).real
         if show:
             _show(vals)
-    np.testing.assert_allclose(vals, ref_vals, rtol=0, atol=1e-7)
+        np.testing.assert_allclose(vals, ref_vals[i], rtol=0, atol=1e-7)
 
 
 def _show(vals):
@@ -105,3 +175,4 @@ def _show(vals):
         print("%.7f, " % v, end="")
         if (i + 1) % 5 == 0:
             print("")
+    print("")

test/harmonic/test_dynamical_matrix.py

@@ -1,10 +1,12 @@
 """Tests for dynamical matrix classes."""
 
+from __future__ import annotations
+
 import numpy as np
 import pytest
 
 from phonopy import Phonopy
-from phonopy.harmonic.dynamical_matrix import DynamicalMatrixGL
+from phonopy.harmonic.dynamical_matrix import DynamicalMatrix, DynamicalMatrixGL
 
 dynmat_ref_000 = [
     0.052897,
@@ -461,7 +463,12 @@ dynmat_ref_555 = [
 @pytest.mark.parametrize(
     "is_compact_fc,lang", [(True, "C"), (False, "C"), (True, "Py"), (False, "Py")]
 )
-def test_dynmat(ph_nacl_nonac, ph_nacl_nonac_compact_fc, is_compact_fc, lang):
+def test_dynmat(
+    ph_nacl_nonac: Phonopy,
+    ph_nacl_nonac_compact_fc: Phonopy,
+    is_compact_fc: bool,
+    lang: str,
+):
     """Test dynamical matrix of NaCl in C and python implementations.
 
     1. Without NAC.
@@ -478,7 +485,7 @@ def test_dynmat(ph_nacl_nonac, ph_nacl_nonac_compact_fc, is_compact_fc, lang):
 
 
 @pytest.mark.parametrize("lang", ["C", "Py"])
-def test_dynmat_dense_svecs(ph_nacl_nonac_dense_svecs: Phonopy, lang):
+def test_dynmat_dense_svecs(ph_nacl_nonac_dense_svecs: Phonopy, lang: str):
     """Test with dense svecs."""
     ph = ph_nacl_nonac_dense_svecs
     dynmat = ph.dynamical_matrix
@@ -563,7 +570,7 @@ def test_dynmat_gonze_lee_short_range_fc(ph_nacl: Phonopy):
     )
 
 
-def test_dynmat_gonze_lee_full_term(ph_nacl):
+def test_dynmat_gonze_lee_full_term(ph_nacl: Phonopy):
     """Test with NAC by Gonze and Lee."""
     dynmat = ph_nacl.dynamical_matrix
     _dynmat = DynamicalMatrixGL(
@@ -576,19 +583,13 @@ def test_dynmat_gonze_lee_full_term(ph_nacl):
     _test_dynmat_252525(_dynmat, dynmat_gonze_lee_full_ref_252525)
 
 
-def test_dynmat_wang(ph_nacl):
+def test_dynmat_wang(ph_nacl_wang: Phonopy):
     """Test with NAC by Wang et al."""
-    nac_params = ph_nacl.nac_params
+    dynmat = ph_nacl_wang.dynamical_matrix
-    nac_params["method"] = "wang"
-    ph_nacl.nac_params = nac_params
-    dynmat = ph_nacl.dynamical_matrix
     _test_dynmat_252525(dynmat, dynmat_wang_ref_252525)
-    # Reset nac_params['method']
-    nac_params.pop("method")
-    ph_nacl.nac_params = nac_params
 
 
-def _test_dynmat(dynmat, lang=None):
+def _test_dynmat(dynmat: DynamicalMatrix, lang=None):
     dtype_complex = "c%d" % (np.dtype("double").itemsize * 2)
     if lang:
         dynmat.run([0, 0, 0], lang=lang)
@@ -609,7 +610,7 @@ def _test_dynmat(dynmat, lang=None):
     np.testing.assert_allclose(dynmat.dynamical_matrix, dynmat_ref, atol=1e-5)
 
 
-def _test_dynmat_252525(dynmat, dynmat_ref, lang=None):
+def _test_dynmat_252525(dynmat: DynamicalMatrix, dynmat_ref: list, lang=None):
     dtype_complex = "c%d" % (np.dtype("double").itemsize * 2)
     if lang:
         dynmat.run([0.25, 0.25, 0.25], lang=lang)
@@ -617,7 +618,9 @@ def _test_dynmat_252525(dynmat, dynmat_ref, lang=None):
         dynmat.run([0.25, 0.25, 0.25])
     # for row in dynmat.dynamical_matrix:
     #     print("".join(["%f, %f, " % (c.real, c.imag) for c in row]))
-    dynmat_ref = (
+
-        np.array(dynmat_ref, dtype="double").view(dtype=dtype_complex).reshape(6, 6)
+    np.testing.assert_allclose(
+        dynmat.dynamical_matrix,
+        np.array(dynmat_ref, dtype="double").view(dtype=dtype_complex).reshape(6, 6),
+        atol=1e-5,
     )
-    np.testing.assert_allclose(dynmat.dynamical_matrix, dynmat_ref, atol=1e-5)

test/phonon/test_group_velocity.py

@@ -39,6 +39,41 @@ def test_gv_nacl(ph_nacl: Phonopy):
     #     print("".join(["%.8f, " % v for v in line]))
 
 
+def test_gv_nacl_wang(ph_nacl_wang: Phonopy):
+    """Test of GroupVelocity by NaCl with Wang's NAC method.
+
+    This test should pass _get_dD_analytical.
+
+    """
+    gv_ref = [
+        14.56800976,
+        14.56800976,
+        14.56800976,
+        14.56800976,
+        14.56800976,
+        14.56800976,
+        25.16351730,
+        25.16351730,
+        25.16351730,
+        1.51378156,
+        1.51378156,
+        1.51378156,
+        1.51378156,
+        1.51378156,
+        1.51378156,
+        -7.84946438,
+        -7.84946438,
+        -7.84946438,
+    ]
+    gv = GroupVelocity(
+        ph_nacl_wang.dynamical_matrix, symmetry=ph_nacl_wang.primitive_symmetry
+    )
+    gv.run([[0.1, 0.1, 0.1]])
+    # for line in gv.group_velocities[0]:
+    #     print("".join(["%.8f, " % v for v in line]))
+    np.testing.assert_allclose(gv.group_velocities[0].ravel(), gv_ref, atol=1e-5)
+
+
 def test_gv_si(ph_si: Phonopy):
     """Test of GroupVelocity by Si.