code: def set_tensor_symmetry_PJ(force_constants,
                             lattice,
                             positions,
                             symmetry):
      
  
      rotations = symmetry.get_symmetry_operations()['rotations']
      translations = symmetry.get_symmetry_operations()['translations']
      symprec = symmetry.get_symmetry_tolerance()
  
      N = len(rotations)
  
      mapa = _get_atom_indices_by_symmetry(lattice,
                                           positions,
                                           rotations,
                                           translations,
                                           symprec)
      cart_rot = np.array([similarity_transformation(lattice, rot).T
                           for rot in rotations])
      cart_rot_inv = np.array([np.linalg.inv(rot) for rot in cart_rot])
      fcm = np.array([force_constants[mapa[n], :, :, :][:, mapa[n], :, :]
                      for n in range(N)])
      s = np.transpose(np.array([np.dot(cart_rot[n],
                                        np.dot(fcm[n], cart_rot_inv[n]))
                                 for n in range(N)]), (0, 2, 3, 1, 4))
      force_constants[:] = np.array(np.average(s, axis=0),
                                    dtype='double',
                                    order='C')
docs: Full force constants are symmetrized using crystal symmetry.
  
      This method extracts symmetrically equivalent sets of atomic pairs and
      take sum of their force constants and average the sum.
  
      Since get_force_constants_disps may include crystal symmetry, this method
      is usually meaningless.