|
| 1 | +--- |
| 2 | +jupytext: |
| 3 | + text_representation: |
| 4 | + extension: .md |
| 5 | + format_name: myst |
| 6 | + format_version: 0.13 |
| 7 | + jupytext_version: 1.10.3 |
| 8 | +kernelspec: |
| 9 | + display_name: Python 3 |
| 10 | + language: python |
| 11 | + name: python3 |
| 12 | +--- |
| 13 | + |
| 14 | +# IFCs Embedding |
| 15 | + |
| 16 | +This section explains how to obtain the phonon modes of a defect system using the IFC (Interatomic Force Constant) embedding approach. This method enables the calculation of both defect-localized and bulk-like phonon modes within a unified framework. We illustrate the process using the Sr[Li$_2$Al$_2$O$_2$N$_2$]:Eu$^{2+}$ example, first computing pristine and defect phonons, then performing the embedding. |
| 17 | + |
| 18 | +```{note} |
| 19 | +The simulation parameters below are for demonstration only and are not converged. |
| 20 | +``` |
| 21 | + |
| 22 | +## 1. Pristine Phonons |
| 23 | + |
| 24 | +The following script creates an AbiPy workflow to compute the phonons of the pristine system using DFPT. It uses the primitive structure of Sr[Li$_2$Al$_2$O$_2$N$_2$] and computes phonons on a 2x2x2 q-mesh. |
| 25 | + |
| 26 | +```python |
| 27 | +import sys |
| 28 | +import os |
| 29 | +import abipy.abilab as abilab |
| 30 | +import abipy.data as abidata |
| 31 | +import abipy.core.structure as structure |
| 32 | +from abipy import flowtk |
| 33 | + |
| 34 | +def make_scf_input(): |
| 35 | + """ |
| 36 | + This function constructs the input file for the GS calculation: |
| 37 | + """ |
| 38 | + pseudodir='pseudos' |
| 39 | + |
| 40 | + pseudos = ('Eu.xml', |
| 41 | + 'Sr.xml', |
| 42 | + 'Al.xml', |
| 43 | + 'Li.xml', |
| 44 | + 'O.xml', |
| 45 | + 'N.xml') |
| 46 | + stru = structure.Structure.from_file("SALON_prim.cif") |
| 47 | + |
| 48 | + # Initialize the input |
| 49 | + gs_inp = abilab.AbinitInput(stru, pseudos=pseudos,pseudo_dir=pseudodir) |
| 50 | + |
| 51 | + # Set variables |
| 52 | + gs_inp.set_vars( |
| 53 | + ecut=10, |
| 54 | + pawecutdg=20, |
| 55 | + diemac=5, |
| 56 | + nstep=100, |
| 57 | + tolvrs=1e-15, |
| 58 | + nband=60, |
| 59 | + nbdbuf=10, |
| 60 | + ) |
| 61 | + |
| 62 | + gs_inp.set_kmesh(ngkpt=[2,2,2],shiftk=[[0.5,0.5,0.5]]) |
| 63 | + |
| 64 | + return gs_inp |
| 65 | + |
| 66 | +def build_flow(options): |
| 67 | + # Working directory (default is the name of the script with '.py' removed and "run_" replaced by "flow_") |
| 68 | + if not options.workdir: |
| 69 | + options.workdir = os.path.basename(sys.argv[0]).replace(".py", "").replace("run_", "flow_") |
| 70 | + flow = flowtk.Flow(options.workdir, manager=options.manager) |
| 71 | + scf_input = make_scf_input() |
| 72 | + phonon_work = flowtk.works.PhononWork.from_scf_input(scf_input,qpoints=[2,2,2],is_ngqpt=True, with_becs=True,ddk_tolerance={"tolwfr":1e-10}) |
| 73 | + flow.register_work(phonon_work) |
| 74 | + return flow |
| 75 | + |
| 76 | +@flowtk.flow_main |
| 77 | +def main(options): |
| 78 | + """ |
| 79 | + This is our main function that will be invoked by the script. |
| 80 | + flow_main is a decorator implementing the command line interface. |
| 81 | + Command line args are stored in `options`. |
| 82 | + """ |
| 83 | + return build_flow(options) |
| 84 | + |
| 85 | +if __name__ == "__main__": |
| 86 | + sys.exit(main()) |
| 87 | +``` |
| 88 | + |
| 89 | +## 2. Defect Phonons |
| 90 | + |
| 91 | +This script computes the phonons of the defect system using finite differences. It uses the supercell defect structure from the $\Delta$SCF calculation and then uses a phonopy workflow. The supercell size is set to [1,1,1], which is equivalent to a $\Gamma$ point calculation. |
| 92 | + |
| 93 | +```python |
| 94 | +import sys |
| 95 | +import os |
| 96 | +import abipy.abilab as abilab |
| 97 | +import abipy.data as abidata |
| 98 | +from abipy.core.structure import Structure |
| 99 | +from abipy import flowtk |
| 100 | +from abipy.flowtk.abiphonopy import PhonopyWork |
| 101 | + |
| 102 | +def make_scf_input(): |
| 103 | + # extract the structure from the deltaSCF calculation. |
| 104 | + structure = Structure.from_file("flow_deltaSCF/w0/t0/outdata/out_GSR.nc") |
| 105 | + pseudodir='pseudos' |
| 106 | + |
| 107 | + pseudos = ('Eu.xml', |
| 108 | + 'Sr.xml', |
| 109 | + 'Al.xml', |
| 110 | + 'Li.xml', |
| 111 | + 'O.xml', |
| 112 | + 'N.xml') |
| 113 | + |
| 114 | + gs_scf_inp = abilab.AbinitInput(structure, pseudos=pseudos,pseudo_dir=pseudodir) |
| 115 | + gs_scf_inp.set_vars(ecut=10, |
| 116 | + pawecutdg=20, |
| 117 | + chksymbreak=0, |
| 118 | + diemac=5, |
| 119 | + prtwf=-1, |
| 120 | + nstep=100, |
| 121 | + toldfe=1e-6, |
| 122 | + chkprim=0, |
| 123 | + ) |
| 124 | + |
| 125 | + # Set DFT+U and spinat parameters according to chemical symbols. |
| 126 | + symb2luj = {"Eu": {"lpawu": 3, "upawu": 7, "jpawu": 0.7}} |
| 127 | + |
| 128 | + gs_scf_inp.set_usepawu(usepawu=1, symb2luj=symb2luj) |
| 129 | + |
| 130 | + n_val = gs_scf_inp.num_valence_electrons |
| 131 | + n_cond = round(20) |
| 132 | + |
| 133 | + spin_up_gs = f"\n{int((n_val - 7) / 2)}*1 7*1 {n_cond}*0" |
| 134 | + spin_up_ex = f"\n{int((n_val - 7) / 2)}*1 6*1 0 1 {n_cond - 1}*0" |
| 135 | + spin_dn = f"\n{int((n_val - 7) / 2)}*1 7*0 {n_cond}*0" |
| 136 | + |
| 137 | + nsppol = 2 |
| 138 | + shiftk = [0, 0, 0] |
| 139 | + ngkpt = [1, 1, 1] |
| 140 | + |
| 141 | + # Build SCF input for the ground state configuration. |
| 142 | + gs_scf_inp.set_kmesh_nband_and_occ(ngkpt, shiftk, nsppol, [spin_up_gs, spin_dn]) |
| 143 | + # Build SCF input for the excited configuration. |
| 144 | + exc_scf_inp = gs_scf_inp.deepcopy() |
| 145 | + exc_scf_inp.set_kmesh_nband_and_occ(ngkpt, shiftk, nsppol, [spin_up_ex, spin_dn]) |
| 146 | + |
| 147 | + return gs_scf_inp |
| 148 | + |
| 149 | + |
| 150 | +def build_flow(options): |
| 151 | + |
| 152 | + if not options.workdir: |
| 153 | + options.workdir = os.path.basename(sys.argv[0]).replace(".py", "").replace("run_", "flow_") |
| 154 | + flow = flowtk.Flow(options.workdir, manager=options.manager) |
| 155 | + |
| 156 | + # Build input for GS calculation |
| 157 | + scf_input = make_scf_input() |
| 158 | + |
| 159 | + # We only need the phonons at Gamma point, so we set the supercell to [1,1,1]. |
| 160 | + work = PhonopyWork.from_gs_input(scf_input, scdims=[1, 1, 1]) |
| 161 | + flow.register_work(work) |
| 162 | + |
| 163 | + return flow |
| 164 | + |
| 165 | + |
| 166 | +@flowtk.flow_main |
| 167 | +def main(options): |
| 168 | + """ |
| 169 | + This is our main function that will be invoked by the script. |
| 170 | + flow_main is a decorator implementing the command line interface. |
| 171 | + Command line args are stored in `options`. |
| 172 | + """ |
| 173 | + return build_flow(options) |
| 174 | + |
| 175 | + |
| 176 | +if __name__ == "__main__": |
| 177 | + sys.exit(main()) |
| 178 | + |
| 179 | +``` |
| 180 | + |
| 181 | +## 3. IFC Embedding: Step-by-Step |
| 182 | + |
| 183 | +The IFC embedding procedure combines the results of the two previous calculations to produce a new `Phonopy` object that contains the phonon modes of the defect system in a large supercell. The embedded phonons are typically created with the `Embedded_phonons.from_phonopy_instances` method : |
| 184 | + |
| 185 | +```{code-cell} |
| 186 | +from abipy.embedding.embedding_ifc import Embedded_phonons |
| 187 | +help(Embedded_phonons.from_phonopy_instances) |
| 188 | +``` |
| 189 | + |
| 190 | + |
| 191 | +### 3.1. Load Pristine and Defect Phonon Data |
| 192 | + |
| 193 | +First, load the DDB file (pristine phonons) and the phonopy object for the defect system: |
| 194 | + |
| 195 | +```{code-cell} |
| 196 | +from abipy.dfpt.converters import ddb_ucell_to_phonopy_supercell |
| 197 | +from abipy.embedding.embedding_ifc import Embedded_phonons |
| 198 | +from pymatgen.io.phonopy import get_pmg_structure |
| 199 | +from abipy.core.kpoints import kmesh_from_mpdivs |
| 200 | +from abipy.abilab import abiopen |
| 201 | +import phonopy |
| 202 | +import numpy as np |
| 203 | +
|
| 204 | +# Load pristine phonons (DFPT, 2x2x2 q-mesh) |
| 205 | +ddb_pristine = abiopen("../workflows_data/flow_phonons/w0/outdata/out_DDB") |
| 206 | +
|
| 207 | +# Load defect phonons (finite difference, supercell) |
| 208 | +ph_defect = phonopy.load( |
| 209 | + supercell_filename="../workflows_data/flow_phonons_doped/w0/outdata/POSCAR", |
| 210 | + force_sets_filename="../workflows_data/flow_phonons_doped/w0/outdata/FORCE_SETS" |
| 211 | +) |
| 212 | +``` |
| 213 | + |
| 214 | +### 3.2. Fold Pristine IFCs to Supercell |
| 215 | + |
| 216 | +The pristine DDB file is first interpolated using `anaget_interpolated_ddb` and then the folding procedure is done with `ddb_ucell_to_phonopy_supercell`: |
| 217 | + |
| 218 | +```{code-cell} |
| 219 | +# Define the large phonon supercell size for the embedding |
| 220 | +sc_size = [2, 2, 4] |
| 221 | +
|
| 222 | +# Generate the list of q-points corresponding to the supercell |
| 223 | +qpts = kmesh_from_mpdivs(mpdivs=sc_size, shifts=[0, 0, 0], order="unit_cell") |
| 224 | +
|
| 225 | +# Interpolate the pristine DDB to obtain force constants on the q-point mesh |
| 226 | +ddb_pristine_inter = ddb_pristine.anaget_interpolated_ddb(qpt_list=qpts) |
| 227 | +
|
| 228 | +# Fold the phonons, and convert to Phonopy object. |
| 229 | +ph_pristine = ddb_ucell_to_phonopy_supercell(ddb_pristine_inter) |
| 230 | +``` |
| 231 | + |
| 232 | +### 3.3. Prepare Structures and Defect Mapping |
| 233 | + |
| 234 | +To embed the IFCs, the algorithm needs to know which atom(s) are substituted or modified. This requires: |
| 235 | +- The defect structure *without* relaxation (i.e., before local relaxation around the defect) |
| 236 | +- The coordinates of the defect atom in both the defect and pristine supercells |
| 237 | + |
| 238 | +```{code-cell} |
| 239 | +# Create the defect structure without relaxation |
| 240 | +structure_defect_wo_relax = ddb_pristine.structure.copy() |
| 241 | +structure_defect_wo_relax.make_supercell([1, 1, 2]) |
| 242 | +structure_defect_wo_relax.replace(0, 'Eu') |
| 243 | +
|
| 244 | +# Index and coordinates of the defect atom (manually identified) |
| 245 | +idefect_defect_stru = 0 |
| 246 | +main_defect_coords_in_defect = structure_defect_wo_relax.cart_coords[idefect_defect_stru] |
| 247 | +
|
| 248 | +idefect_pristine_stru = 0 # (manually identified) |
| 249 | +main_defect_coords_in_pristine = get_pmg_structure(ph_pristine.supercell).cart_coords[idefect_pristine_stru] |
| 250 | +``` |
| 251 | + |
| 252 | +### 3.4. Run the Embedding Algorithm |
| 253 | + |
| 254 | +The `Embedded_phonons.from_phonopy_instances` method combines the pristine and defect phonon data. The algorithm: |
| 255 | +- Maps atoms between the pristine and defect supercells (crucial step) |
| 256 | +- Extracts IFCs from both systems |
| 257 | +- Replaces pristine IFCs with defect IFCs within a cutoff radius around the defect |
| 258 | +- Retains pristine IFCs elsewhere |
| 259 | +- Enforces the Acoustic Sum Rule (ASR) |
| 260 | + |
| 261 | +You can print details of the modifications by setting `verbose=True`. |
| 262 | + |
| 263 | +```{code-cell} |
| 264 | +emb_ph = Embedded_phonons.from_phonopy_instances( |
| 265 | + phonopy_pristine=ph_pristine, |
| 266 | + phonopy_defect=ph_defect, |
| 267 | + structure_defect_wo_relax=structure_defect_wo_relax, |
| 268 | + main_defect_coords_in_defect=main_defect_coords_in_defect, |
| 269 | + main_defect_coords_in_pristine=main_defect_coords_in_pristine, |
| 270 | + substitutions_list=[[idefect_pristine_stru, "Eu"]], |
| 271 | + cut_off_mode="auto" |
| 272 | +) |
| 273 | +``` |
| 274 | + |
| 275 | +If the structure mapping fails (e.g., due to mismatched coordinates), the code will notify you: |
| 276 | + |
| 277 | +```{code-cell} |
| 278 | +emb_ph_failed = Embedded_phonons.from_phonopy_instances( |
| 279 | + phonopy_pristine=ph_pristine, |
| 280 | + phonopy_defect=ph_defect, |
| 281 | + structure_defect_wo_relax=structure_defect_wo_relax, |
| 282 | + main_defect_coords_in_defect=main_defect_coords_in_defect, |
| 283 | + main_defect_coords_in_pristine=main_defect_coords_in_pristine + np.array([0.1, 0.0, 0.0]), |
| 284 | + substitutions_list=[[idefect_pristine_stru, "Eu"]], |
| 285 | + cut_off_mode="auto", |
| 286 | + verbose=False |
| 287 | +) |
| 288 | +``` |
| 289 | + |
| 290 | +Finally, note that the `Embedded_phonons` class provides methods for further analysis and data export. You can use `get_gamma_freq_with_vec_abipy_fmt()` to compute the $\Gamma$-point phonon frequencies and eigenvectors in AbiPy format, or `to_ddb()` to convert the embedded phonons back to an Abinit DDB file. Also, you will find in `abipy.embedding.utils_ifc` function to computethe localization ratio of the phonon modes as well as helper function to draw phonon eigenvectors with VESTA. |
| 291 | + |
| 292 | + |
| 293 | + |
| 294 | +```{note} |
| 295 | +For additional examples of the use of this module, especially for the embedding of different defect types (vacancy, interstitial), you can check the tests examples located in `abipy/embedding/tests/test_embedding_ifc.py`. |
| 296 | +``` |
| 297 | + |
| 298 | + |
| 299 | + |
| 300 | + |
| 301 | + |
| 302 | + |
| 303 | + |
0 commit comments