nac_term = (A1 (x) A2) / B * coef.

def _nac(self, q_direction):
        """nac_term = (A1 (x) A2) / B * coef.
        """
        num_atom = self._pcell.get_number_of_atoms()
        nac_q = np.zeros((num_atom, num_atom, 3, 3), dtype='double')
        if (np.abs(q_direction) < 1e-5).all():
            return nac_q

        rec_lat = np.linalg.inv(self._pcell.get_cell())
        nac_factor = self._dynmat.get_nac_factor()
        Z = self._dynmat.get_born_effective_charges()
        e = self._dynmat.get_dielectric_constant()
        q = np.dot(rec_lat, q_direction)

        B = self._B(e, q)
        for i in range(num_atom):
            A_i = self._A(q, Z, i)
            for j in range(num_atom):
                A_j = self._A(q, Z, j)
                nac_q[i, j] = np.outer(A_i, A_j) / B

        num_satom = self._scell.get_number_of_atoms()
        N = num_satom // num_atom

        return nac_q * nac_factor / N