code: def _parse_fc(self, f, natom, dim):
          
  
          ndim = np.prod(dim)
          fc = np.zeros((natom, natom * ndim, 3, 3), dtype='double', order='C')
          for k, l, i, j in np.ndindex((3, 3, natom, natom)):
              line = f.readline()
              for i_dim in range(ndim):
                  line = f.readline()
                  # fc[i, j * ndim + i_dim, k, l] = float(line.split()[3])
                  fc[j, i * ndim + i_dim, l, k] = float(line.split()[3])
          return fc
docs: Parse force constants part
  
          Physical unit of force cosntants in the file is Ry/au^2.