Bugfix on symmetries in the supercell for systems with 48 symmetries …

…in the point group

cellconstructor/Phonons.py

@@ -13,7 +13,7 @@
 import os, sys
 import scipy, scipy.optimize
 
-
+import itertools
 import cellconstructor.Structure as Structure
 import cellconstructor.symmetries as symmetries
 import cellconstructor.Methods as Methods
@@ -1438,6 +1438,64 @@ def GetRamanVector(self, pol_in, pol_out):
 
         return v
 
+    def GetRamanActive(self, use_spglib = False):
+        """
+        This simple subroutines tries to guess by symmetry analisys which mode is active or not.
+        If a raman tensor is present, it will be used to test the activity, otherwise, a random one will
+        be generated
+
+        Parameters
+        ----------
+            use_spblib: bool
+                If True the spglib library is used to initialize symmetries.
+                Usefull if the phonon matrix is in a super cell.
+
+        Results
+        -------
+            raman_activity_mask : ndarray(size = (3*nat), dtype = bool)
+                A mask that is False or True if a mode in the unit cell is Raman-active or not.
+        """
+
+        there_is_raman_tensor = True
+        if self.raman_tensor is None:
+            there_is_raman_tensor = False 
+            self.raman_tensor = np.zeros((3,3, 3*self.structure.N_atoms), dtype = np.double)
+            self.raman_tensor[:,:,:] = np.random.uniform( size = self.raman_tensor.shape)
+
+            # Get the symmetries
+            qe_sym = symmetries.QE_Symmetry(self.structure)
+            if use_spglib:
+                qe_sym.SetupFromSPGLIB()
+            else:
+                qe_sym.SetupQPoint()
+
+            # Symmetrize the effective charges
+            qe_sym.ApplySymmetryToRamanTensor(self.raman_tensor)
+
+        # Save a debugging one
+        self.save_qe("Raman")
+
+        # Simulate the Raman signal for all possible incoming and outcoming polarizations
+        res = np.zeros(self.structure.N_atoms * 3, dtype = np.double)
+        for i, j in itertools.product(range(3) , range(3)):
+            pol_in = np.zeros(3)
+            pol_in[i] = 1
+            pol_out = np.zeros(3)
+            pol_out[j] = 1
+
+            res += self.GetRamanResponce(pol_in, pol_out)
+
+        print("total_raman_res:", res)
+
+        is_raman_active = res > 1e-5    
+
+        # Delete the random raman tensor if any
+        if not there_is_raman_tensor:
+            self.raman_tensor = None
+
+        return is_raman_active
+
+
 
     def GenerateSupercellDyn(self, supercell_size):
         """
@@ -2341,8 +2399,14 @@ def SymmetrizeSupercell(self, supercell_size = None):
 
         for iq, q in enumerate(self.q_tot):
             self.dynmats[iq] = fcq[iq, :, :]
+
+        # Symmetrize also the effective charges and the Raman Tensor if any
+        if not self.effective_charges is None:
+            qe_sym.ApplySymmetryToEffCharge(self.eff_charges)
+        if not self.raman_tensor is None:
+            qe_sym.ApplySymmetryToRamanTensor(self.raman_tensor)
 
-    def Symmetrize(self, verbose = False, asr = "custom"):
+    def Symmetrize(self, verbose = False, asr = "custom", use_spglib = False):
         """
         SYMMETRIZE THE DYNAMICAL MATRIX
         ===============================
@@ -2357,7 +2421,12 @@ def Symmetrize(self, verbose = False, asr = "custom"):
             asr : string
                 The kind of the acustic sum rule. Allowed are 'crystal', 'simple' or 'custom'.
                 for crystal and simple refer to the quantum-espresso guide.
+            use_spglib : bool
+                If True, the simmetrization is performed with SPGLIB in the supercell
         """
+
+        if use_spglib:
+            self.SymmetrizeSupercell()
 
         qe_sym = symmetries.QE_Symmetry(self.structure)
 
@@ -2366,6 +2435,13 @@ def Symmetrize(self, verbose = False, asr = "custom"):
 
         for iq,q in enumerate(self.q_tot):
             self.dynmats[iq] = fcq[iq, :, :]
+
+        # Symmetrize also the effective charges and the Raman Tensor if any
+        if not self.effective_charges is None:
+            qe_sym.ApplySymmetryToEffCharge(self.effective_charges)
+        if not self.raman_tensor is None:
+            qe_sym.ApplySymmetryToRamanTensor(self.raman_tensor)
+
 
 
     def ApplySumRule(self, kind = "custom"):
@@ -2451,7 +2527,7 @@ def GetIRActive(self, use_spglib = False):
 
         # Simulate the IR signal
         Ir = self.GetIRIntensities()
-        is_ir_active = Ir > 1e-9    
+        is_ir_active = Ir > 1e-8    
 
         # Delete the random effective charges if added
         if not there_are_eff_charges:

cellconstructor/Structure.py

@@ -951,6 +951,28 @@ def save_bcs(self, filename, symmetry_file = ""): # STILL NOT WORKING
 
         fp.close()
 
+    def get_xcoords(self):
+        """
+        Returns the crystalline coordinates
+        """
+
+        assert self.has_unit_cell
+
+        xcoords = np.zeros(self.coords)
+        for i in range(self.N_atoms):
+            xcoords[i,:] = Methods.covariant_coordinates(self.unit_cell, self.coords[i,:])
+
+        return xcoords
+    def set_from_xcoords(self, xcoords):
+        """
+        Set the cartesian coordinates from crystalline
+        """
+
+        assert self.has_unit_cell
+
+        for i in range(self.N_atoms):
+            self.coords[i,:]  = self.unit_cell.T.dot(xcoords[i,:])
+
 
     def save_scf(self, filename, alat = 1, avoid_header=False):
         """
@@ -1024,6 +1046,17 @@ def fix_coords_in_unit_cell(self):
 
         # Delete duplicate atoms
         self.delete_copies()
+
+    def fix_wigner_seitz(self):
+        """
+        Atoms will be replaced in the periodic images inside the wigner_seitz cell
+        """
+
+        assert self.has_unit_cell, "Error, the wigner_seitz is defined for periodic boundary conditions"
+
+        for i in range(self.N_atoms):
+            new_r = Methods.get_closest_vector(self.unit_cell, self.coords[i,:])
+            self.coords[i, :] = new_r
 
     def get_strct_conventional_cell(self):
         """
@@ -1143,7 +1176,12 @@ def get_itau(self, unit_cell_structure):
         ----------
             - unit_cell_structure : Structure()
                 The structure of the unit cell used to generate this supercell structure.
-                
+
+        Results
+        -------
+            - itau : ndarray (size = nat_sc, type = int)
+                For each atom in the supercell contains the index of the corrisponding
+                atom in the unit_cell, starting from 1 to unit_cell_structure.N_atoms (included)
         """
 
         itau = np.zeros( self.N_atoms, dtype = np.intc)
@@ -1159,6 +1197,15 @@ def get_itau(self, unit_cell_structure):
 
         return itau
 
+    def get_sublattice_vectors(self, unit_cell_structure):
+        """
+        Get the lattice vectors that connects the atom of this supercell structure to those of
+        the unit_cell structure.
+        """
+
+        itau = self.get_itau(unit_cell_structure) - 1 
+        return self.coords[:,:] - unit_cell_structure.coords[itau[:], :]
+
     def generate_supercell(self, dim, itau = None, QE_convention = True, get_itau = False):
         """
         This method generate a supercell of specified dimension, replicating the system

cellconstructor/symmetries.py

@@ -357,6 +357,84 @@ def ApplySymmetryToEffCharge(self, eff_charges):
         for i in range(nat):
             eff_charges[i, :, :] = Methods.convert_matrix_cart_cryst(new_eff_charges[i, :, :], self.QE_at.T, True)
 
+    def ApplySymmetryToRamanTensor(self, raman_tensor):
+        """
+        SYMMETRIZE RAMAN TENSOR
+        ============================
+
+        This subroutine applies the symmetries to the raman tensor
+
+        As always, the raman_tensor will be modified by this subroutine.
+
+        Parameters
+        ----------
+            - raman_tensor : ndarray (size = (3, 3, 3*nat))
+                The raman tensor. The first two indices indicate
+                the polarization of the incoming/outcoming field, while the last one
+                is the atomic/cartesian coordinate
+        """
+
+        pol1, pol2, at_cart = np.shape(raman_tensor)
+
+        assert pol1 == pol2 
+        assert pol2 == 3
+        assert at_cart == 3*self.QE_nat, "Error, the structure and effective charges are not compatible"
+
+        # Apply the permutation on the electric fields
+        raman_tensor += np.einsum("abc->bac", raman_tensor)
+        raman_tensor /= 2
+
+        # Apply the sum rule
+        # The sum over all the atom for each cartesian coordinate should be zero.
+        rt_reshaped = raman_tensor.reshape((3,3,self.QE_nat, 3))
+
+        # Sum over all the atomic indices
+        tot_sum = np.sum(rt_reshaped, axis = 2)
+
+        # Rebuild the shift to the tensor of the correct shape
+        shift = np.tile(tot_sum, (self.QE_nat, 1, 1, 1))
+
+        # Place the number of atoms at the correct position
+        # From the first to the third
+        shift = np.einsum("abcd->bcad", shift)
+
+        # Now we apply the sum rule
+        rt_reshaped -= shift / self.QE_nat
+        new_tensor = np.zeros(np.shape(rt_reshaped), dtype = np.double)
+
+        # Get the raman tensor in crystal components
+        for i in range(self.QE_nat):
+            rt_reshaped[:,:, i, :] = Methods.convert_3tensor_to_cryst(rt_reshaped[:,:, i, :], self.QE_at.T)
+
+        # Apply translations
+        if self.QE_translation_nr > 1:
+            for i in range(self.QE_translation_nr):
+                irt = self.QE_translations_irt[:, i] - 1
+                for j in range(self.QE_nat):
+                    new_mat = rt_reshaped[:,:, irt[j], :]
+                    new_tensor += new_mat
+
+            rt_reshaped = new_tensor / self.QE_translation_nr
+            new_tensor[:,:,:,:] = 0.
+
+        # Apply rotations
+        for i in range(self.QE_nsym):
+            irt = self.QE_irt[i, :] - 1
+
+            for j in range(self.QE_nat):
+                # Apply the symmetry to the 3 order tensor
+                new_mat = np.einsum("ai, bj, ck, ijk", self.QE_s[:,:,i], self.QE_s[:,:,i], self.QE_s[:,:,i], rt_reshaped[:,:, irt[j], :])
+                #new_mat = self.QE_s[:,:, i].dot( eff_charges[irt[j], :, :].dot(self.QE_s[:,:,i].T))
+                new_tensor[:,:,j,:] += new_mat
+
+        new_tensor /= self.QE_nsym
+
+        # Convert back into cartesian
+        for i in range(self.QE_nat):
+            rt_reshaped[:, :, i, :] = Methods.convert_3tensor_to_cryst(new_tensor[:,:,i,:], self.QE_at.T, True)
+
+        # Compress again the notation
+        raman_tensor[:,:,:] = rt_reshaped.reshape((3,3, 3*self.QE_nat))
 
 
     def ApplySymmetryToSecondOrderEffCharge(self, dM_drdr, apply_asr = True):
@@ -996,6 +1074,12 @@ def SetupFromSPGLIB(self):
             # Extract the rotation and the fractional translation
             rot = sym[:,:3]
 
+            # Check if the rotation is equal to the first one
+            if np.sum( (rot - symmetries[0][:,:3])**2 ) < 0.1 and n_syms == 0 and i > 0:
+                # We got all the rotations
+                n_syms = i 
+                break
+
             # Extract the point group
             if n_syms == 0:
                 self.QE_s[:,:, i] = rot.T
@@ -1004,15 +1088,10 @@ def SetupFromSPGLIB(self):
                 irt = GetIRT(self.structure, sym)
                 self.QE_irt[i, :] = irt + 1 #Py to Fort
 
-            # Check if the rotation is equal to the first one
-            if np.sum( (rot - symmetries[0][:,:3])**2 ) < 0.1 and n_syms == 0 and i > 0:
-                # We got all the rotations
-                n_syms = i 
-                break
 
         if n_syms == 0:
             n_syms = len(symmetries)
-
+        
         # From the point group symmetries, get the supercell
         n_supercell = len(symmetries) // n_syms
         self.QE_translation_nr = n_supercell

requirements.txt

@@ -1,3 +1,4 @@
+setuptools
 numpy
 scipy
 ase

requirements2.txt

@@ -0,0 +1,5 @@
+setuptools
+numpy
+scipy
+ase==3.16.0
+spglib

scripts/cellconstructor_test.py

@@ -242,6 +242,23 @@ def setUp(self):
         self.dynSnSe = DownloadDynSnSe()
         self.dynSky = DownloadDynSky()
 
+        # Test a simple NaCl setup
+        a = 5.41
+        self.NaCl = CC.Structure.Structure(2)
+        self.NaCl.atoms = ["Na","Cl"]
+        self.NaCl.coords[1,:] = (a/2, a/2, a/2)
+        self.NaCl.unit_cell[0,:] = (a/2, a/2, 0)
+        self.NaCl.unit_cell[1,:] = (a/2, 0, a/2)
+        self.NaCl.unit_cell[2,:] = (0, a/2, a/2)
+        self.NaCl.has_unit_cell = True
+        self.NaCl.set_masses({"Na": 22.989769, "Cl": 35.453})
+        self.NaCldyn = CC.Phonons.Phonons(self.NaCl)
+        self.NaCldyn.dynmats[0][:,:] = np.random.uniform(size = (6, 6))
+        self.NaCldyn.dynmats[0] += self.NaCldyn.dynmats[0].T
+        self.NaCldyn.Symmetrize()
+
+
+
 
 
     def test_generate_supercell(self):
@@ -281,6 +298,39 @@ def test_IR_modes(self):
         """
         raise NotImplementedError("Error, this test must be implemented")
 
+    def test_IR_Raman_on_NaCl(self):
+        """
+        We generate a simple NaCl structure and predict using symmetries if the mode is IR or Raman active
+        """
+        self.setUp()
+        self.NaCldyn.Symmetrize()
+
+        for i in range(100):
+            get_ir_activity = self.NaCldyn.GetIRActive()
+            assert (get_ir_activity == [False, False, False, True, True, True]).all(), get_ir_activity
+
+        for i in range(100):
+            get_raman_activity = self.NaCldyn.GetRamanActive()
+            assert not get_raman_activity.any()
+
+        # Check if the activity vector is actually good
+        self.NaCldyn.dynmats[0][:,:] = np.random.uniform(size = (3*self.NaCl.N_atoms, 3*self.NaCl.N_atoms))
+        self.NaCldyn.effective_charges = np.random.uniform(size = (self.NaCl.N_atoms, 3, 3))
+        self.NaCldyn.raman_tensor = np.random.uniform(size = (3, 3, 3 * self.NaCl.N_atoms))
+
+        self.NaCldyn.Symmetrize()
+
+        # Check that the raman tensor is really zero
+        assert np.max(np.abs(self.NaCldyn.raman_tensor)) < 1e-10
+
+        # Get the IR responce
+        ir_responce = self.NaCldyn.GetIRIntensities()
+
+        # Check that if the intensities along the degenerate modes are all equal
+        assert np.max(np.abs(ir_responce[3:] - ir_responce[3])) < 1e-10
+
+
+
 
     def test_change_phonon_cell(self):
         """