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.

def set_tensor_symmetry_PJ(force_constants,
                           lattice,
                           positions,
                           symmetry):
    """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.

    """

    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')