Fixes + new plotting functions (#308)

* fix asr, + add the nac params of interstitial. * combine plotting fonctions for the 1D case and multi-D case. + plot eigenenergies of the four states in delta SCF. * fix asr, + add the nac params of interstitial. * combine plotting fonctions for the 1D case and multi-D case. + plot eigenenergies of the four states in delta SCF. * Solve compatibility issue with new phonopy compact_fc_to_full_fc() function. + Adapt the ifc embedding tests with the correction of the asr formula * Small improvements for PR --------- Co-authored-by: Julien Bouquiaux <[email protected]>
abinit · Feb 4, 2025 · a4a6c59 · a4a6c59
1 parent 9cd89c5
commit a4a6c59
Show file tree

Hide file tree

Showing 8 changed files with 462 additions and 209 deletions.
diff --git a/abipy/dfpt/converters.py b/abipy/dfpt/converters.py
@@ -771,11 +771,17 @@ def ddb_ucell_to_phonopy_supercell(unit_ddb=None,unit_ddb_filepath=None,nac=True
     phonopy_ucell= unit_ddb.anaget_phonopy_ifc(ngqpt=ngqpt, asr=1, dipdip=0, chneut=1,
                                           set_masses=True)
 
-    # fc from (uc_size x sc_size) to (sc_size x sc_size)
 
-    full_fc = force_constants.compact_fc_to_full_fc(phonon=phonopy_ucell,
+    # fc from (uc_size x sc_size) to (sc_size x sc_size)
+    try:
+        full_fc = force_constants.compact_fc_to_full_fc(primitive=phonopy_ucell.primitive,
+                                                    compact_fc=phonopy_ucell.force_constants)
+
+    except TypeError: #old compact_fc_to_full_fc function (phonopy <= 2.32)
+        full_fc = force_constants.compact_fc_to_full_fc(phonon=phonopy_ucell,  
                                                     compact_fc=phonopy_ucell.force_constants)
 
+
     # create a phonopy object with supercell structure
     phonopy_supercell = Phonopy(unitcell=phonopy_ucell.supercell,# the new unit cell is the 'old' supercell
                                 supercell_matrix=[1, 1, 1],  # sup_size= unit_size, gamma only

diff --git a/abipy/embedding/embedding_ifc.py b/abipy/embedding/embedding_ifc.py
@@ -197,15 +197,15 @@ def from_phonopy_instances(cls,
                         print(f"Diff IFC = \n {ifc_pristine[i][j]-ifc_defect[mapping.index(i)][mapping.index(j)]}" )
 
                     ifc_emb[i][j]=factor_ifc*ifc_defect[mapping.index(i)][mapping.index(j)]
-
-                # enforce ASR
-        if asr==True:
+                    
+        # enforce ASR
+        if asr:
             print(f"\n Enforce ASR")
-            for alpha in [0,1,2]:
-                for beta in [0,1,2]:
-                    for i,atom1 in enumerate(stru_emb):
-                        ifc_emb[i][i][alpha,beta] = -(accoustic_sum(ifc_emb,i)[alpha,beta]-ifc_emb[i][i][alpha,beta])
-
+            sum_ac=np.sum(ifc_emb,axis=1)
+            for i,atom1 in enumerate(stru_emb):
+                for alpha in [0,1,2]:
+                    ifc_emb[i][i][alpha][alpha] = - (sum_ac[i][alpha][alpha]-ifc_emb[i][i][alpha][alpha])
+        
         ########################
         # change the nac params
         ########################
@@ -216,6 +216,11 @@ def from_phonopy_instances(cls,
 
             if vacancies_list is not None:
                 nac_params_emb["born"]=np.delete(nac_params_emb["born"],vacancies_list,0)
+            if interstitial_list is not None:
+                for interstial in interstitial_list:
+                    nac=np.append(nac_params_emb["born"],(stru_emb[-1].specie.common_oxidation_states[0])*np.eye(3))
+
+                    nac_params_emb["born"]=nac.reshape(len(stru_emb),3,3)
         else:
             nac_params_emb=None
 

diff --git a/abipy/embedding/tests/test_embedding_ifc.py b/abipy/embedding/tests/test_embedding_ifc.py
@@ -13,16 +13,11 @@
 from abipy.embedding.utils_ifc import localization_ratio
 from abipy.embedding.embedding_ifc import Embedded_phonons
 from pymatgen.io.phonopy import get_pmg_structure,get_phonopy_structure
-
-
-try:
-    import phonopy
-except ImportError:
-    Phonopy = None
+import phonopy
 
 
 class Embedding_ifcTest(AbipyTest):
-
+    
     def test_embedding_vacancy(self):
 
         self.skip_if_not_phonopy()
@@ -31,12 +26,16 @@ def test_embedding_vacancy(self):
         # open a 4x4x4 q-mesh ddb of CaO
         ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/CaO_444_DDB"))
 
-        # interpolate on a 3x3x3 q-mesh
+        ## interpolate on a 3x3x3 q-mesh
+
         qpts=kmesh_from_mpdivs(mpdivs=[3,3,3],shifts=[0,0,0],order="unit_cell")
         ddb_pristine_333=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts)
 
+
+
         ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_pristine_333,nac=True)
 
+
         ddb_defect = DdbFile(abidata.ref_file("refs/embedding_ifc/CaO_16at_vacancy_DDB"))
         ph_defect=ddb_defect.anaget_phonopy_ifc()
 
@@ -54,7 +53,7 @@ def test_embedding_vacancy(self):
         structure_defect_wo_relax.remove_sites(indices=[idefect_defect_stru])
         structure_defect_wo_relax.sort()
 
-        #
+        # 
         idefect_pristine_stru=27
         main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell).cart_coords[idefect_pristine_stru]
 
@@ -64,29 +63,31 @@ def test_embedding_vacancy(self):
                                phonopy_defect=ph_defect,
                                structure_defect_wo_relax=structure_defect_wo_relax,
                                main_defect_coords_in_pristine=main_defect_coords_in_pristine,
-                               main_defect_coords_in_defect=main_defect_coords_in_defect,
-                               substitutions_list=None,
-                               vacancies_list=[27],
-                               interstitial_list=None,
+                               main_defect_coords_in_defect=main_defect_coords_in_defect, 
+                               substitutions_list=None, 
+                               vacancies_list=[27],       
+                               interstitial_list=None,  
                                factor_ifc=3, #  increasing ifc, this will induce a fictious local mode above the bulk frequencies
                                cut_off_mode='auto',
                                verbose=0,
                                asr=True,)
-
+        
         #test the conversion to ddb
-        ph_emb.to_ddb(tmp_dir + "out_DDB")
-        ddb_emb = DdbFile(tmp_dir + "out_DDB")
+        ph_emb.to_ddb(tmp_dir+"out_DDB")
+        ddb_emb = DdbFile(tmp_dir+"out_DDB")
 
         modes=ddb_emb.anaget_phmodes_at_qpoint([0,0,0])
         freqs_anaddb=modes.phfreqs[0]
         freqs_phonopy,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt()
 
-        self.assert_almost_equal(freqs_anaddb[-1],7.31410150e-02,decimal=5)
-        self.assert_almost_equal(freqs_phonopy[-1],7.31411352e-02,decimal=5)
-        self.assert_almost_equal(max(abs(freqs_anaddb-freqs_phonopy)),0,decimal=5)
+        self.assert_almost_equal(freqs_anaddb[-1],0.07314101496974906,decimal=5)
+        self.assert_almost_equal(freqs_phonopy[-1],0.07182725446338956,decimal=5)
+        self.assert_almost_equal(max(abs(freqs_anaddb-freqs_phonopy)),0.0020214752349544465,decimal=5)
+
+        ratio=localization_ratio(vecs)
+
+        self.assert_almost_equal(ratio[-1],7.129688982510354,decimal=5)
 
-        ratio = localization_ratio(vecs)
-        self.assert_almost_equal(ratio[-1], 7.41692860, decimal=5)
 
     def test_embedding_substitution(self):
 
@@ -96,20 +97,24 @@ def test_embedding_substitution(self):
         ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/SrCl2_DDB"))
 
         ## interpolate on a 6x6x6 q-mesh
+
         qpts=kmesh_from_mpdivs(mpdivs=[6,6,6],shifts=[0,0,0],order="unit_cell")
         ddb_pristine_666=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts)
 
+
+
         ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_pristine_666,nac=False)
 
-        # test with phonopy loading
+        # test with phonopy loading 
         ph_defect = phonopy.load(supercell_filename=abidata.ref_file("refs/embedding_ifc/SrCl2_Eu_POSCAR"),
                          force_sets_filename=abidata.ref_file("refs/embedding_ifc/SrCl2_Eu_FORCE_SETS"))
 
+
         ########
         # We need first to create the defect structure without relax
         structure_defect_wo_relax=ddb_pristine.structure.copy()
         structure_defect_wo_relax.make_supercell(3)
-        structure_defect_wo_relax.replace(0, 'Eu')
+        structure_defect_wo_relax.replace(0,'Eu')
         structure_defect_wo_relax.sort()
 
         # index of the sub. = 26 (in defect structure), this is found manually
@@ -120,58 +125,68 @@ def test_embedding_substitution(self):
         idefect_pristine_stru=0
         main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell).cart_coords[idefect_pristine_stru]
 
+
         ph_emb=Embedded_phonons.from_phonopy_instances(
                                phonopy_pristine=ph_pristine,
                                phonopy_defect=ph_defect,
                                structure_defect_wo_relax=structure_defect_wo_relax,
                                main_defect_coords_in_pristine=main_defect_coords_in_pristine,
-                               main_defect_coords_in_defect=main_defect_coords_in_defect,
-                               substitutions_list=[[idefect_pristine_stru,"Eu"]],
-                               vacancies_list=None,
-                               interstitial_list=None,
+                               main_defect_coords_in_defect=main_defect_coords_in_defect, 
+                               substitutions_list=[[idefect_pristine_stru,"Eu"]], 
+                               vacancies_list=None,       
+                               interstitial_list=None,  
                                cut_off_mode='auto',
                                verbose=0,
                                asr=True,)
 
-        freqs, vecs = ph_emb.get_gamma_freq_with_vec_abipy_fmt()
-        self.assert_almost_equal(freqs[-1], 0.026440349386134484, decimal=5)
 
-        ratio = localization_ratio(vecs)
-        self.assert_almost_equal(ratio[446], 21.65624596304, decimal=5)
+        freqs,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt()
+
+        self.assert_almost_equal(freqs[-1],0.026431844680473826,decimal=5)
+
+        ratio=localization_ratio(vecs)
+
+        self.assert_almost_equal(ratio[446],2.354695997610233,decimal=5)
+
 
     def test_embedding_interstitial(self):
 
         self.skip_if_not_phonopy()
         tmp_dir = tempfile.mkdtemp()
+
         ddb_pristine = DdbFile(abidata.ref_file("refs/embedding_ifc/C_conv_DDB"))
 
+
         qpts=kmesh_from_mpdivs(mpdivs=[3,3,3],shifts=[0,0,0],order="unit_cell")
         ddb_unit_333=ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts,)
 
-        ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_unit_333, nac=False,)
+        ph_pristine=ddb_ucell_to_phonopy_supercell(unit_ddb=ddb_unit_333,
+                                                   nac=False,)
 
-        # test with phonopy loading
+
+        # test with phonopy loading 
         ddb_defect = DdbFile(abidata.ref_file("refs/embedding_ifc/C_interstitial_DDB"))
         ph_defect=ddb_defect.anaget_phonopy_ifc()
 
         ########
+
         ### model of the split-interstitial as one vacancy + 2 inter.
 
-        stru= abiopen(abidata.ref_file("refs/embedding_ifc/C_sc.cif"))
+        stru=abiopen(abidata.ref_file("refs/embedding_ifc/C_sc.cif"))
         stru.remove_sites(indices=[14])
-        stru.append(species="C", coords=[0.600, 0.5000, 0.2500])
-        stru.append(species="C", coords=[0.400, 0.5000, 0.2500])
-        structure_defect_wo_relax = stru
-
-        # main defect is
+        stru.append(species="C",coords=[0.600, 0.5000, 0.2500])
+        stru.append(species="C",coords=[0.400, 0.5000, 0.2500])
+        structure_defect_wo_relax=stru
+        
+        # main defect is 
         main_defect_coords_in_pristine=get_pmg_structure(ph_pristine.supercell)[31].coords
         main_defect_coords_in_defect=main_defect_coords_in_pristine
 
         ph_emb=Embedded_phonons.from_phonopy_instances(phonopy_pristine=ph_pristine,
                                                phonopy_defect=ph_defect,
                                                structure_defect_wo_relax=structure_defect_wo_relax,
                                                main_defect_coords_in_pristine=main_defect_coords_in_pristine,
-                                               main_defect_coords_in_defect=main_defect_coords_in_defect,
+                                               main_defect_coords_in_defect=main_defect_coords_in_defect, 
                                                vacancies_list=[31],
                                                interstitial_list=[['C',[4.2885, 3.5737, 1.7869]],
                                                                   ['C',[2.8590, 3.5737, 1.7869]]
@@ -182,9 +197,12 @@ def test_embedding_interstitial(self):
                                                asr=True)
 
 
+        freqs,vecs=ph_emb.get_gamma_freq_with_vec_abipy_fmt()
+
+        self.assert_almost_equal(freqs[-1],0.21212027299739988,decimal=5)
+
+        ratio=localization_ratio(vecs)
+
+        self.assert_almost_equal(ratio[650],78.29179108928511,decimal=5)
 
-        freqs, vecs = ph_emb.get_gamma_freq_with_vec_abipy_fmt()
-        self.assert_almost_equal(freqs[-1], 0.212826358519, decimal=5)
 
-        ratio = localization_ratio(vecs)
-        self.assert_almost_equal(ratio[650], 78.76445591, decimal=5)