code: def _nac(self, q_direction):
          
          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
docs: nac_term = (A1 (x) A2) / B * coef.