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_old(force_constants,
                            lattice,  # column vectors
                            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()

    fc_bak = force_constants.copy()

    # Create mapping table between an atom and the symmetry operated atom
    # map[ i, j ]
    # i: atom index
    # j: operation index
    mapping = []
    for pos_i in positions:
        map_local = []
        for rot, trans in zip(rotations, translations):
            rot_pos = np.dot(pos_i, rot.T) + trans
            for j, pos_j in enumerate(positions):
                diff = pos_j - rot_pos
                diff -= np.rint(diff)
                diff = np.dot(diff, lattice.T)
                if np.linalg.norm(diff) < symprec:
                    map_local.append(j)
                    break
        mapping.append(map_local)
    mapping = np.array(mapping)

    # Look for the symmetrically equivalent force constant tensors
    for i, pos_i in enumerate(positions):
        for j, pos_j in enumerate(positions):
            tmp_fc = np.zeros((3, 3), dtype='double')
            for k, rot in enumerate(rotations):
                cart_rot = similarity_transformation(lattice, rot)

                # Reverse rotation of force constant is summed
                tmp_fc += similarity_transformation(cart_rot.T,
                                                    fc_bak[mapping[i, k],
                                                           mapping[j, k]])
            # Take average and set to new force cosntants
            force_constants[i, j] = tmp_fc / len(rotations)