GenerateSupercell for the dynamical matrix fix

cellconstructor/Phonons.py

@@ -1885,7 +1885,7 @@ def GetRamanActive(self, use_spglib = False):
 
 
 
-    def GenerateSupercellDyn(self, supercell_size, img_thr = 1e-6):
+    def GenerateSupercellDyn(self, supercell_size=None, img_thr = 1e-6):
         """
         GENERATE SUPERCEL DYN
         =====================
@@ -1899,14 +1899,18 @@ def GenerateSupercellDyn(self, supercell_size, img_thr = 1e-6):
         ----------
             supercell_size : array int (size=3)
                 the dimension of the cell on which you want to generate the new
-                Phonon
+                Phonon. If None, it is inferred from the q points.
 
         Results
         -------
             dyn_supercell : Phonons()
                 A Phonons class of the supercell
 
         """
+        if supercell_size is None:
+            supercell_size = self.GetSupercell()
+
+
         # First check if the q vectors are compatible with the supercell
         if not symmetries.CheckSupercellQ(self.structure.unit_cell, supercell_size, self.q_tot):
             print("Q points:", self.q_tot)
@@ -1916,10 +1920,31 @@ def GenerateSupercellDyn(self, supercell_size, img_thr = 1e-6):
 
         super_struct = self.structure.generate_supercell(supercell_size)
 
+        nat = self.structure.N_atoms
+        nat_sc = super_struct.N_atoms
+
         dyn_supercell = Phonons(super_struct, nqirr = 1, force_real = True)
 
         dyn_supercell.dynmats[0] = self.GetRealSpaceFC(supercell_size, img_thr = img_thr)
 
+        # Check also effective charges and dielectric tensor
+        if self.dielectric_tensor is not None:
+            dyn_supercell.dielectric_tensor = self.dielectric_tensor
+        if self.effective_charges is not None:
+            itau = super_struct.get_itau(self.structure) - 1
+            dyn_supercell.effective_charges = np.zeros((nat_sc, 3, 3), dtype = np.double)
+            for i in range(nat_sc):
+                i_uc = itau[i]
+                dyn_supercell.effective_charges[i, :, :] = self.effective_charges[i_uc, :, :]
+        if self.raman_tensor is not None:
+            itau = super_struct.get_itau(self.structure) - 1
+            dyn_supercell.raman_tensor = np.zeros((3, 3, 3*nat_sc), dtype = np.double)
+            for i in range(3 * nat_sc):
+                i_at = i // 3
+                j_coord = i % 3
+                i_uc = itau[i_at]
+                dyn_supercell.raman_tensor[:, :, i] = self.raman_tensor[:, :, i_uc + j_coord]
+
         return dyn_supercell
 
 

cellconstructor/Structure.py

@@ -1349,6 +1349,9 @@ def get_ase_atoms(self):
         """
         This method returns the ase atoms structure, ready for computations.
 
+        This function automatically transform isotopes like deuterium 'D' in 'H' atoms,
+        otherwise ase crashes.
+
         Results
         -------
             - atoms : ase.Atoms()
@@ -1362,7 +1365,10 @@ def get_ase_atoms(self):
         # Get thee atom list
         atm_list = []
         for i in range(self.N_atoms):
-            atm_list.append(ase.Atom(self.atoms[i], self.coords[i,:]))
+            atm_lbl = self.atoms[i]
+            if atm_lbl == "D":
+                atm_lbl = "H"
+            atm_list.append(ase.Atom(atm_lbl, self.coords[i,:]))
 
         atm = ase.Atoms(atm_list)
 
@@ -2329,4 +2335,4 @@ def get_structures_from_phonopy_supercells(ph_supercells):
 
         all_structures.append(s)
 
-    return all_structures
+    return all_structures

setup.py

@@ -72,7 +72,7 @@
 
 
 setup( name = "CellConstructor",
-       version = "1.4.1",
+       version = "1.5.0",
        description = "Python utilities that is interfaced with ASE for atomic crystal analysis",
        author = "Lorenzo Monacelli",
        url = "https://github.com/mesonepigreco/CellConstructor",

tests/TestEffChargeInterp/test_effective_charges_supercell.py

@@ -0,0 +1,42 @@
+import cellconstructor as CC
+import cellconstructor.Phonons 
+import cellconstructor.ForceTensor
+
+
+import sys, os
+import time
+import numpy as np 
+import pytest
+
+def test_effective_charges_supercell():
+    total_path = os.path.dirname(os.path.abspath(__file__))
+    os.chdir(total_path)
+
+    # Load the dynamical matrix
+    dyn = CC.Phonons.Phonons("dyn", 32)
+
+    # Get enforce the ASR on the effective charges
+    dyn.Symmetrize()
+    ef_charge = dyn.effective_charges.copy()
+
+    super_dyn = dyn.GenerateSupercellDyn()
+
+    # Check the effective charges ASR 
+    for i in range(3):
+        asr_check = np.sum(ef_charge[:, i, :])
+
+        assert np.abs(asr_check) < 1e-6, "ASR not enforced for the uc effective charges on component %d" % i
+
+        # Check also the new effective charges
+        asr_check = np.sum(super_dyn.effective_charges[:, i, :])
+        assert np.abs(asr_check) < 1e-6, "ASR not enforced for the super cell effective charges on component %d" % i
+
+
+    # Check the effective charges
+    assert np.allclose(ef_charge, super_dyn.effective_charges[:dyn.structure.N_atoms, :, :]), "Effective charges not correctly copied"
+
+    print("Correctly copied the effective charges")
+
+
+if __name__ == "__main__":
+    test_effective_charges_supercell()