Include Jochym's new implementation of force constants symmetrization

README.md

@@ -2,3 +2,4 @@ phonopy
 =======
 
 Phonon code
+For the details, see http://phonopy.sourceforge.net/

doc/changelog.rst

@@ -3,13 +3,13 @@
 Change Log
 ==========
 
-Feb-10-2014: Version 1.9.3
+Feb-10-2015: Version 1.9.3
 -----------------------------
 
 * ``MAGMOM`` tag is now available (:ref:`magmom_tag`).
 * Spglib update.
 
-Jan-4-2014: Version 1.9.2
+Jan-4-2015: Version 1.9.2
 -----------------------------
 
 * Behaviors of ``--wien2k``, ``--abinit``, ``--pwscf`` options have

phonopy/harmonic/force_constants.py

@@ -269,6 +269,38 @@ def set_tensor_symmetry(force_constants,
     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.
+    """
+    N = len(rotations)
+    K = len(positions)
+    
+    rpos = np.dot(positions, np.transpose(rotations, (0, 2, 1))) + translations
+    diff = positions - np.transpose(np.tile(rpos, (K, 1, 1, 1)), (2, 1, 0, 3))
+    m = (abs(diff - diff.round()) < symprec).all(axis=-1)
+    mapa = np.tile(np.arange(K, dtype='intc'), (K,1))
+    mapa = np.array([mapa[mr] for mr in m])
+    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')
+
+def set_tensor_symmetry_old(force_constants,
+                        lattice,
+                        positions,
+                        rotations,
+                        translations,
+                        symprec):
+    """
+    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.
     """