diff --git a/cellconstructor/Phonons.py b/cellconstructor/Phonons.py index f099d13c..4ddf8e6c 100644 --- a/cellconstructor/Phonons.py +++ b/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,7 +1899,7 @@ 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 ------- @@ -1907,6 +1907,10 @@ def GenerateSupercellDyn(self, supercell_size, img_thr = 1e-6): 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 diff --git a/cellconstructor/Structure.py b/cellconstructor/Structure.py index c6539039..bcb1f916 100644 --- a/cellconstructor/Structure.py +++ b/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 \ No newline at end of file + return all_structures diff --git a/setup.py b/setup.py index 26dfb917..646087dd 100644 --- a/setup.py +++ b/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", diff --git a/tests/TestEffChargeInterp/test_effective_charges_supercell.py b/tests/TestEffChargeInterp/test_effective_charges_supercell.py new file mode 100644 index 00000000..793a8971 --- /dev/null +++ b/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()