Skip to content

Commit

Permalink
Solve compatibility issue with new phonopy compact_fc_to_full_fc() fu…
Browse files Browse the repository at this point in the history 
…nction.

+ Adapt the ifc embedding tests with the correction of the asr formula
  • Loading branch information
Julien Bouquiaux committed Jan 29, 2025
1 parent 88c2a3e commit 484704f
Show file tree
Hide file tree
Showing 3 changed files with 71 additions and 42 deletions.
10 changes: 8 additions & 2 deletions abipy/dfpt/converters.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
101 changes: 62 additions & 39 deletions abipy/embedding/tests/test_embedding_ifc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,13 +16,13 @@


try:
import phonopy
import phonopy
except ImportError:
Phonopy = None


class Embedding_ifcTest(AbipyTest):

def test_embedding_vacancy(self):

self.skip_if_not_phonopy()
Expand All @@ -31,12 +31,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()

Expand All @@ -54,7 +58,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]

Expand All @@ -64,29 +68,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):

Expand All @@ -96,20 +102,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
Expand All @@ -120,58 +130,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]]
Expand All @@ -182,9 +202,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)
2 changes: 1 addition & 1 deletion requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ plotly
ipython
chart-studio
click
phonopy<=2.31.2
phonopy
seekpath
ase
# TODO remove after https://github.com/materialsproject/emmet/issues/768 is fixed
Expand Down

0 comments on commit 484704f

Please sign in to comment.