1""" 2Converters between abinit/abipy format and other external tools. 3Some portions of the code have been imported from the ConvertDDB.py script 4developed by Hu Xe, Eric Bousquet and Aldo Romero. 5""" 6 7import os 8import itertools 9import warnings 10import numpy as np 11from monty.dev import requires 12 13from pymatgen.io.phonopy import get_phonopy_structure, get_pmg_structure 14from pymatgen.io.vasp.inputs import Poscar 15import abipy.core.abinit_units as abu 16from abipy.dfpt.ddb import DdbFile 17from abipy.abio.factories import minimal_scf_input 18from abipy.electrons.gsr import GsrFile 19from monty.os import makedirs_p 20try: 21 from phonopy import Phonopy 22 from phonopy.file_IO import write_FORCE_CONSTANTS, parse_FORCE_CONSTANTS, parse_BORN, parse_FORCE_SETS 23 #from phonopy.interface.phonopy_yaml import PhonopyYaml 24 from phonopy.interface.calculator import get_default_physical_units, get_force_constant_conversion_factor 25except ImportError: 26 Phonopy = None 27 28 29@requires(Phonopy, "phonopy not installed!") 30def abinit_to_phonopy(anaddbnc, supercell_matrix, symmetrize_tensors=False, output_dir_path=None, 31 prefix_outfiles="", symprec=1e-5, set_masses=False): 32 """ 33 Converts the interatomic force constants(IFC), born effective charges(BEC) and dielectric 34 tensor obtained from anaddb to the phonopy format. Optionally writes the 35 standard phonopy files to a selected directory: FORCE_CONSTANTS, BORN (if BECs are available) 36 POSCAR of the unit cell, POSCAR of the supercell. 37 38 The conversion is performed taking the IFC in the Wigner–Seitz supercell with weights 39 as produced by anaddb and reorganizes them in a standard supercell multiple of the 40 unit cell. Operations are vectorized using numpy. This may lead to the allocation of 41 large arrays in case of very large supercells. 42 43 Performs a check to verify if the two codes identify the same symmetries and it gives a 44 warning in case of failure. Mismatching symmetries may lead to incorrect conversions. 45 46 Args: 47 anaddbnc: an instance of AnaddbNcFile. Should contain the output of the IFC analysis, 48 the BEC and the dielectric tensor. 49 supercell_matrix: the supercell matrix used for phonopy. Any choice is acceptable, however 50 the best agreement between the abinit and phonopy results is obtained if this is set to 51 a diagonal matrix with on the diagonal the ngqpt used to generate the anaddb.nc. 52 symmetrize_tensors: if True the tensors will be symmetrized in the Phonopy object and 53 in the output files. This will apply to IFC, BEC and dielectric tensor. 54 output_dir_path: a path to a directory where the phonopy files will be created 55 prefix_outfiles: a string that will be added as a prefix to the name of the written files 56 symprec: distance tolerance in Cartesian coordinates to find crystal symmetry in phonopy. 57 It might be that the value should be tuned so that it leads to the the same symmetries 58 as in the abinit calculation. 59 set_masses: if True the atomic masses used by abinit will be added to the PhonopyAtoms 60 and will be present in the returned Phonopy object. This should improve compatibility 61 among abinit and phonopy results if frequencies needs to be calculated. 62 63 Returns: 64 An instance of a Phonopy object that contains the IFC, BEC and dieletric tensor data. 65 """ 66 67 ifc = anaddbnc.ifc 68 nac_params = None 69 becs = None 70 epsinf = None 71 if anaddbnc.becs is not None and anaddbnc.epsinf is not None: 72 becs = anaddbnc.becs.values 73 epsinf = anaddbnc.epsinf 74 75 # according to the phonopy website 14.399652 is not the coefficient for abinit 76 # probably it relies on the other conventions in the output. 77 nac_params = {"born": becs, "dielectric": epsinf, "factor": 14.399652} 78 79 s = anaddbnc.structure 80 81 phon_at = get_phonopy_structure(s) 82 if set_masses: 83 phon_at.masses = [anaddbnc.amu[n] for n in phon_at.numbers] 84 85 # use phonopy to get the proper supercell given by the primitive and the matrix 86 # and convert it to pymatgen 87 phonon = Phonopy(phon_at, supercell_matrix, primitive_matrix=np.eye(3), nac_params=nac_params, 88 symprec=symprec) 89 phon_supercell = phonon.get_supercell() 90 supercell = get_pmg_structure(phon_supercell) 91 92 abi_hall_num = s.abi_spacegroup.get_spglib_hall_number() 93 spglib_hall_num = phonon.symmetry.dataset["hall_number"] 94 if abi_hall_num != spglib_hall_num: 95 warnings.warn("The hall number obtained based on the DDB symmetries differs " 96 f"from the one calculated with spglib: {abi_hall_num} versus " 97 f"{spglib_hall_num}. The conversion may be incorrect. Try changing symprec.") 98 99 # convert to phonopy units 100 at_cart = ifc.atoms_cart_coord * abu.Bohr_Ang 101 ifccc = ifc.ifc_cart_coord * abu.Ha_eV / abu.Bohr_Ang ** 2 102 weights = ifc.ifc_weights 103 latt = supercell.lattice 104 105 ifcph = np.zeros((len(s), len(supercell), 3, 3)) 106 107 # loop over the atoms in the primitive cell 108 # other operations are vectorized using numpy arrays. Some array may require large allocations 109 for i, (site, c_list, w_list) in enumerate(zip(s, at_cart, weights)): 110 111 ind_w = np.where(w_list > 0) 112 ifccc_loc = ifccc[i, ind_w[0]] 113 114 w_list = w_list[ind_w] 115 c_list = c_list[ind_w] 116 117 # align the coordinates of the first atom in the list (the site under consideration) 118 # with the site in the primitive cell. 119 c_list = c_list - c_list[0] + site.coords 120 121 # convert to fractional coordinates as needed by the Lattice to get the distances 122 f_list = latt.get_fractional_coords(c_list) 123 sc_fcoords = supercell.frac_coords 124 125 # construct the list of sites of the supercell that are closer to sites in 126 # the primitive cell 127 dist_and_img = [latt.get_distance_and_image(f_list[0], fc) for fc in sc_fcoords] 128 # the function gives the translation of the image, but it should be applied to the coordinates. 129 # Only the positions are needed 130 nearest_sc_fcoords = [fc + trasl for (_, trasl), fc in zip(dist_and_img, sc_fcoords)] 131 132 # divide by the corresponding weights. Elements with weights 0 were discarded above 133 ifccc_loc = np.transpose(ifccc_loc, (0, 2, 1)) / w_list[:, None, None] 134 135 # create an array with all the possible pairs 136 # instantiating this array seems slow but seems still faster than the required loops 137 coord_pairs = np.array(list(itertools.product(nearest_sc_fcoords, f_list))) 138 139 # find the pairs that match between the coordinates of the modified supercell and the f_list 140 ind_match = np.where(np.abs(coord_pairs[:, 0] - coord_pairs[:, 1]).sum(axis=1) < 1e-6)[0] 141 # set the ifc for phonopy in the final array corresponding to the matching indices. 142 n_points_f_list = len(f_list) 143 ifcph[i, ind_match // n_points_f_list] = ifccc_loc[ind_match % n_points_f_list] 144 145 phonon.set_force_constants(ifcph) 146 if symmetrize_tensors: 147 phonon.symmetrize_force_constants() 148 149 if output_dir_path: 150 makedirs_p(output_dir_path) 151 152 fc_filepath = os.path.join(output_dir_path, prefix_outfiles+"FORCE_CONSTANTS") 153 write_FORCE_CONSTANTS(phonon.get_force_constants(), fc_filepath) 154 155 if becs is not None and epsinf is not None: 156 born_filepath = os.path.join(output_dir_path, prefix_outfiles+"BORN") 157 write_BORN(phon_at, borns=becs, epsilon=epsinf, filename=born_filepath, 158 symmetrize_tensors=symmetrize_tensors) 159 160 poscar_filepath = os.path.join(output_dir_path, prefix_outfiles+"POSCAR") 161 poscar = Poscar(s) 162 poscar.write_file(poscar_filepath, significant_figures=15) 163 164 supercell_filepath = os.path.join(output_dir_path, prefix_outfiles+"supercell_POSCAR") 165 superce_poscar = Poscar(supercell) 166 superce_poscar.write_file(supercell_filepath, significant_figures=15) 167 168 return phonon 169 170 171@requires(Phonopy, "phonopy not installed!") 172def phonopy_to_abinit(unit_cell, supercell_matrix, out_ddb_path, ngqpt=None, qpt_list=None, 173 force_constants=None, force_sets=None, born=None, 174 primitive_matrix="auto", symprec=1e-5, tolsym=None, supercell=None, 175 calculator=None, manager=None, workdir=None, pseudos=None, verbose=False): 176 """ 177 Converts the data from phonopy to an abinit DDB file. The data can be provided 178 in form of arrays or paths to the phonopy files that should be parsed. 179 The minimal input should contains the FORCE_CONSTANTS or FORCE_SETS. 180 If BORN is present the Born effective charges (BEC) and dielectric 181 tensor will also be added to the DDB. 182 183 The best agreement is obtained with supercell_matrix and ngqpt being 184 equivalent (i.e. supercell_matrix a diagonal matrix with ngqpt as diagonal 185 elements). Non diagonal supercell_matrix are allowed as well, but the information 186 encoded in the DDB will be the result of an interpolation done through phonopy. 187 188 Phonopy is used to convert the IFC to the dynamical matrix. However, in order to 189 determine the list of q-points in the irreducible Brillouin zone and to prepare the 190 base for the final DDB file, abinit will be called for a very short and inexpensive run. 191 192 Performs a check to verify if the two codes identify the same symmetries and it gives a 193 warning in case of failure. Mismatching symmetries may lead to incorrect conversions. 194 195 Args: 196 unit_cell: a |Structure| object that identifies the unit cell used for the phonopy 197 calculation. 198 supercell_matrix: a 3x3 array representing the supercell matrix used to generated the 199 forces with phonopy. 200 out_ddb_path: a full path to the file where the new DDB will be written 201 ngqpt: a list of 3 elements indicating the grid of q points that will be used in the DDB. 202 qpt_list: alternatively to ngqpt an explicit list of q-points can be provided here. 203 At least one among ngqpt and qpt_list should be defined. 204 force_constants: an array with shape (num atoms unit cell, num atoms supercell, 3, 3) 205 containing the force constants. Alternatively a string with the path to the 206 FORCE_CONSTANTS file. This or force_set should be defined. If both given this 207 has precedence. 208 force_sets: a dictionary obtained from the force sets generated with phonopy. 209 Alternatively a string with the path to the FORCE_SETS file. This or force_constants 210 should be defined. 211 born: a dictionary with "dielectric" and "born" keywords as obtained from the nac_params 212 in phonopy. Alternatively a string with the path to the BORN file. Notice that 213 the "factor" attribute is not taken into account, so the values should be in 214 default phonopy units. 215 primitive_matrix: a 3x3 array with the primitive matrix passed to Phonopy. "auto" will 216 use spglib to try to determine it automatically. If the DDB file should contain the 217 actual unit cell this should be the identity matrix. 218 symprec: distance tolerance in Cartesian coordinates to find crystal symmetry in phonopy. 219 It might be that the value should be tuned so that it leads to the the same symmetries 220 as in the abinit calculation. 221 tolsym: Gives the tolerance to identify symmetries in abinit. See abinit documentation for 222 more details. 223 supercell: if given it should represent the supercell used to get the force constants, 224 without any perturbation. It will be used to match it to the phonopy supercell 225 and sort the IFC in the correct order. 226 calculator: a string with the name of the calculator. Will be used to set the conversion 227 factor for the force constants coming from phonopy. 228 manager: |TaskManager| object. If None, the object is initialized from the configuration file 229 pseudos: List of filenames or list of |Pseudo| objects or |PseudoTable| object. It will be 230 used by abinit to generate the base DDB file. If None the abipy.data.hgh_pseudos.HGH_TABLE 231 table will be used. 232 verbose: verbosity level. Set it to a value > 0 to get more information 233 workdir: path to the directory where the abinit calculation will be executed. 234 235 Returns: 236 a DdbFile instance of the file written in out_ddb_path. 237 """ 238 239 if ngqpt is None and qpt_list is None: 240 raise ValueError("at least one among nqgpt and qpt_list should be defined") 241 242 if force_sets is None and force_constants is None: 243 raise ValueError("at least one of force_sets and force_constants should be provided") 244 245 phon_at = get_phonopy_structure(unit_cell) 246 247 if isinstance(force_constants, str): 248 force_constants = parse_FORCE_CONSTANTS(filename=force_constants) 249 elif force_constants is not None: 250 force_constants = np.array(force_constants) 251 force_sets = None 252 253 if isinstance(force_sets, str): 254 force_sets = parse_FORCE_SETS(filename=force_sets) 255 256 # no nac_params here, otherwise they will be used for the interpolation 257 phonon = Phonopy(phon_at, supercell_matrix, primitive_matrix=primitive_matrix, nac_params=None, 258 symprec=symprec, calculator=calculator) 259 260 primitive = get_pmg_structure(phonon.primitive) 261 262 if isinstance(born, str): 263 born = parse_BORN(phonon.primitive, filename=born) 264 265 if supercell is not None: 266 ph_supercell = get_pmg_structure(phonon.supercell) 267 if not np.allclose(supercell.lattice.matrix, ph_supercell.lattice.matrix): 268 raise RuntimeError("The lattice of the supercells do not match") 269 sc_mapping = [] 270 for i, site_orig in enumerate(supercell): 271 for j, site_ph in enumerate(ph_supercell): 272 d = supercell.lattice.get_distance_and_image(site_orig.frac_coords, site_ph.frac_coords)[0] 273 if d < 1e-5: 274 sc_mapping.append(j) 275 break 276 else: 277 raise RuntimeError(f"Could not find a match for site {i} with coords " 278 f"{site_orig.cart_coords} in the supercell.") 279 280 # cross check that the same atom was not matched twice 281 n_matches = len(set(sc_mapping)) 282 if n_matches < len(supercell): 283 raise RuntimeError(f"Found matches for {n_matches} different atoms in the supercell: {sc_mapping}") 284 285 force_constants = force_constants[:, sc_mapping] 286 287 if force_constants is not None: 288 phonon.set_force_constants(force_constants) 289 else: 290 phonon.dataset = force_sets 291 phonon.produce_force_constants() 292 293 if calculator: 294 units = get_default_physical_units(calculator) 295 fc_factor = get_force_constant_conversion_factor(units["force_constants_unit"], None) 296 phonon.set_force_constants(phonon.force_constants * fc_factor) 297 298 if pseudos is None: 299 from abipy.data.hgh_pseudos import HGH_TABLE 300 pseudos = HGH_TABLE 301 302 inp = minimal_scf_input(primitive, pseudos) 303 304 # get the qpoints list if not defined 305 if qpt_list is None: 306 inp["ngkpt"] = ngqpt 307 qpt_list = inp.abiget_ibz(verbose=verbose)[0] 308 309 dm_list = get_dm(phonon, qpt_list, primitive) 310 311 if born is not None: 312 # for the conversion of the BEC the zion (i.e. the ionic charge of the pseudo) 313 # it is an additive factor and should be the same that goes in the header of the DDB, 314 # so take it from the pseudos used to generate it. 315 zion = inp.valence_electrons_per_atom 316 born_data = generate_born_deriv(born, zion, primitive) 317 else: 318 born_data = None 319 320 inp = minimal_scf_input(primitive, pseudos) 321 if tolsym is not None: 322 inp["tolsym"] = tolsym 323 task = inp.run_in_shell(workdir=workdir, manager=manager, verbose=verbose) 324 325 # use the output of abinit to check that the spacegroup identified by 326 # phonopy and abinit are the same. 327 with GsrFile(task.opath_from_ext("GSR.nc")) as gsr: 328 abi_spg = gsr.structure.abi_spacegroup.spgid 329 spglib_spg = phonon.symmetry.dataset["number"] 330 if abi_spg != spglib_spg: 331 warnings.warn("The space group number obtained based on the DDB symmetries differs " 332 f"from the one calculated with spglib: {abi_spg} versus " 333 f"{spglib_spg}. The convertion may be incorrect. Try changing symprec or tolsym.") 334 335 tmp_ddb_path = task.opath_from_ext("DDB") 336 337 ddb = DdbFile(tmp_ddb_path) 338 # remove the blocks generated by the calculation and that are meaningless 339 ddb.remove_block(dord=0) 340 ddb.remove_block(dord=1) 341 342 add_data_ddb(ddb, dm_list, qpt_list, born_data) 343 344 ddb.write(out_ddb_path) 345 346 new_ddb = DdbFile(out_ddb_path) 347 return new_ddb 348 349 350def generate_born_deriv(born, zion, structure): 351 """ 352 Helper function to generate the portion of the derivatives in the DDB 353 that are related to the Born effective charges and dielectric tensor, 354 starting from the data available in phonopy format. 355 356 Args: 357 born: a dictionary with "dielectric" and "born" keywords as obtained from the nac_params 358 in phonopy. 359 zion: the ionic charge of each atom in the system. It should be in the same order 360 as the one present in the header of the DDB. 361 structure: a pymatgen |Structure| of the unit cell. 362 363 Returns: 364 a complex numpy array with shape (len(structure)+2, 3, len(structure)+2, 3). Only the 365 parts relative to the BECs and dielectric tensors will be filled. 366 """ 367 natoms = len(structure) 368 mpert = natoms + 2 # only these perturbations are needed here 369 born_data = np.zeros((3, mpert, 3, mpert), dtype=complex) 370 371 eps_e = born["dielectric"] 372 bec = np.array(born["born"]).transpose((0, 2, 1)) 373 374 # Assume that for the generated DDB acell = [1,1,1] and rprimd = rprim. Should be in Bohr. 375 rprim = structure.lattice.matrix * abu.Ang_Bohr 376 rprimd = rprim 377 378 volume_bohr = np.dot(rprimd[:, 0], np.cross(rprimd[:, 1], rprimd[:, 2])) 379 gprimd = np.linalg.inv(rprimd) 380 dij = np.identity(3) 381 # BEC 382 for ipert1 in range(natoms): # ipert1 is atom position deriv 383 ipert2 = natoms + 1 # E field deriv 384 dm1 = np.matmul(rprimd, 385 np.matmul(bec[ipert1, :, :] - dij[:, :] * zion[ipert1], gprimd)) 386 for idir1 in range(3): 387 for idir2 in range(3): 388 born_data[idir1, ipert1, idir2, ipert2] = dm1[idir1, idir2] * 2 * np.pi + 0.0j 389 born_data[idir2, ipert2, idir1, ipert1] = dm1[idir1, idir2] * 2 * np.pi + 0.0j 390 # epsinf 391 ipert1 = natoms + 1 392 ipert2 = natoms + 1 393 dm1 = np.matmul(gprimd.transpose(), 394 np.matmul(dij[:, :] - eps_e[:, :], gprimd)) 395 for idir1 in range(3): 396 for idir2 in range(3): 397 born_data[idir1, ipert1, idir2, ipert2] = dm1[idir1, idir2] * np.pi * volume_bohr + 0.0j 398 born_data = born_data.transpose((1, 0, 3, 2)) 399 return born_data 400 401 402@requires(Phonopy, "phonopy not installed!") 403def get_dm(phonon, qpt_list, structure): 404 """ 405 Helper function to generate the dynamical matrix in the abinit conventions 406 for a list of q-points from a Phonopy object. 407 408 Args: 409 phonon: a Phonopy object with force constants. 410 qpt_list: a list of fractional coordinates of q-points for which the 411 dynamical matrix should be generated. 412 structure: a pymatgen |Structure| of the primitive cell. 413 414 Returns: 415 a list of arrays with the dynamical matrices of the selected q-points. 416 """ 417 natom = len(structure) 418 rprim = structure.lattice.matrix * abu.Ang_Bohr 419 # assume acell is [1., 1., 1.] 420 rprimd = rprim 421 masses = phonon.masses 422 dm_list = [] 423 mass_tile = np.tile(masses, (3, 1)).T 424 mass_matrix = np.sqrt(np.einsum("ij,kl", mass_tile, mass_tile)) 425 # get the difference in coordinates for each pair of atoms 426 # diff_coords[i,j] is equivalent to structure[i].coords - structure[j].coords 427 coord_matrix = np.repeat(structure.cart_coords[np.newaxis, :, :], natom, axis=0) 428 diff_coords = np.transpose(coord_matrix, (1, 0, 2)) - coord_matrix 429 rlatt = structure.lattice.reciprocal_lattice 430 for q in qpt_list: 431 q_cart = rlatt.get_cartesian_coords(q) 432 # the phase exp(-i * q.r) 433 phase = np.exp(-1j * np.einsum("ijk,k", diff_coords, q_cart)) 434 dm = phonon.get_dynamical_matrix_at_q(q).reshape(natom, 3, natom, 3) 435 # the following is rprim * dm[ipert1,:,ipert2,:] * rprim.T 436 dm = np.einsum("ij,kjlm,nm->kiln", rprimd, dm, rprimd) 437 dm = dm * mass_matrix / phase[:, None, :, None] 438 dm *= abu.eV_Ha / abu.Ang_Bohr**2 439 440 dm_list.append(dm) 441 442 return dm_list 443 444 445def add_data_ddb(ddb, dm_list, qpt_list, born_data): 446 """ 447 Helper function to add the blocks for the dynamical matrix and BECs to a DdbFile object. 448 449 Args: 450 ddb: a DdbFile object to be modified. 451 dm_list: the list of dynamical matrices to be added. 452 qpt_list: the list of q-points corresponding to dm_list. 453 born_data: the data corresponding to BECs and dielectric tensor. If None 454 these part will not be set in the DDB. 455 """ 456 dm_data = {} 457 natom = len(ddb.structure) 458 for q, dm in zip(qpt_list, dm_list): 459 q_data = {} 460 for ipert1 in range(natom): 461 for idir1 in range(3): 462 for ipert2 in range(natom): 463 for idir2 in range(3): 464 q_data[(idir1+1, ipert1+1, idir2+1, ipert2+1)] = dm[ipert1, idir1, ipert2, idir2] 465 466 # for gamma set also the born data if present 467 if np.allclose(q, (0, 0, 0)) and born_data is not None: 468 ipert2 = natom + 1 469 for ipert1 in range(natom): 470 for idir1 in range(3): 471 for idir2 in range(3): 472 q_data[(idir1+1, ipert1+1, idir2+1, ipert2+1)] = born_data[ipert1, idir1, ipert2, idir2] 473 q_data[(idir2+1, ipert2+1, idir1+1, ipert1+1)] = born_data[ipert2, idir2, ipert1, idir1] 474 475 for idir1 in range(3): 476 for idir2 in range(3): 477 q_data[(idir1+1, ipert2+1, idir2+1, ipert2+1)] = born_data[ipert2, idir1, ipert2, idir2] 478 479 dm_data[tuple(q)] = q_data 480 481 ddb.set_2nd_ord_data(dm_data, replace=True) 482 483 484@requires(Phonopy, "phonopy not installed!") 485def tdep_to_abinit(unit_cell, fc_path, supercell_matrix, supercell, out_ddb_path, ngqpt=None, 486 qpt_list=None, primitive_matrix="auto", lotosplitting_path=None, symprec=1e-5, 487 tolsym=None, manager=None, workdir=None, pseudos=None, verbose=False): 488 """ 489 Converts the files produced by TDEP to an abinit DDB file. If the lotosplitting 490 file is provided the BEC and dielectric tensor will also be added to the DDB. 491 492 The conversion is performed by first extracting the force constants in phonopy format 493 and then using phonopy_to_abinit_py to generate the DDB file. See the phonopy_to_abinit_py 494 docstring for further details about the second step of the conversion. 495 Notice that the supercell used by TDEP should be provided in order to properly 496 sort the IFC to match the order required by phonopy. 497 498 Args: 499 unit_cell: a |Structure| object that identifies the unit cell used for the TDEP 500 calculation. 501 fc_path: the path to the forceconstants file produced by TDEP. 502 supercell_matrix: a 3x3 array representing the supercell matrix used to generated the 503 forces with TDEP. 504 supercell: the supercell used by TDEP to get the force constants, without any 505 perturbation (usually named with extension ssposcar). It will be used to match it to 506 the phonopy supercell and sort the IFC in the correct order. 507 out_ddb_path: a full path to the file where the new DDB will be written 508 ngqpt: a list of 3 elements indicating the grid of q points that will be used in the DDB. 509 qpt_list: alternatively to ngqpt an explicit list of q-points can be provided here. 510 At least one among ngqpt and qpt_list should be defined. 511 primitive_matrix: a 3x3 array with the primitive matrix passed to Phonopy. "auto" will 512 use spglib to try to determine it automatically. If the DDB file should contain the 513 actual unit cell this should be the identity matrix. 514 lotosplitting_path: path to the lotosplitting file produced by TDEP. If None no BEC 515 contribution will be set to the DDB. 516 symprec: distance tolerance in Cartesian coordinates to find crystal symmetry in phonopy. 517 It might be that the value should be tuned so that it leads to the the same symmetries 518 as in the abinit calculation. 519 tolsym: Gives the tolerance to identify symmetries in abinit. See abinit documentation for 520 more details. 521 manager: |TaskManager| object. If None, the object is initialized from the configuration file 522 pseudos: List of filenames or list of |Pseudo| objects or |PseudoTable| object. It will be 523 used by abinit to generate the base DDB file. If None the abipy.data.hgh_pseudos.HGH_TABLE 524 table will be used. 525 verbose: verbosity level. Set it to a value > 0 to get more information 526 workdir: path to the directory where the abinit calculation will be executed. 527 528 Returns: 529 a DdbFile instance of the file written in out_ddb_path. 530 """ 531 532 fc = parse_tdep_fc(fc_path, unit_cell, supercell) 533 534 born = None 535 if lotosplitting_path: 536 eps, becs = parse_tdep_lotosplitting(lotosplitting_path) 537 born = {"dielectric": eps, "born": becs, "factor": 1} 538 539 ddb = phonopy_to_abinit(unit_cell=unit_cell, force_constants=fc, supercell_matrix=supercell_matrix, ngqpt=ngqpt, 540 qpt_list=qpt_list, out_ddb_path=out_ddb_path, born=born, pseudos=pseudos, 541 primitive_matrix=primitive_matrix, supercell=supercell, manager=manager, 542 workdir=workdir, symprec=symprec, verbose=verbose, tolsym=tolsym) 543 544 return ddb 545 546 547def parse_tdep_fc(fc_path, unit_cell, supercell): 548 """ 549 Parses a forceconstants file produced by TDEP an converts it to an array in the 550 phonopy format. 551 552 Args: 553 fc_path: path to the forceconstants file 554 unit_cell: a |Structure| object with the unit cell used for the calculation 555 in TDEP. 556 supercell: the supercell used for the calculation in TDEP. 557 558 Returns: 559 a comple numpy array with shape (len(unit_cell), len(supercell), 3, 3) 560 """ 561 natoms = len(unit_cell) 562 fc = np.zeros((natoms, len(supercell), 3, 3)) 563 # parse the fc file and put it in the format required for phonopy 564 # this will be sorted according to the tdep supercell atoms order, not the phonopy supercell. 565 with open(fc_path, "rt") as f: 566 lines = f.readlines() 567 u_latt = unit_cell.lattice 568 sc_latt = supercell.lattice 569 iline = 2 570 for iat in range(natoms): 571 n_neighbours = int(lines[iline].split()[0]) 572 iline += 1 573 for jn in range(n_neighbours): 574 iprim = int(lines[iline].split()[0]) - 1 575 r = [float(sp) for sp in lines[iline + 1].split()] 576 # find matching atom in the supercell 577 fcoords = sc_latt.get_fractional_coords(u_latt.get_cartesian_coords(r) + unit_cell[iprim].coords) 578 for isc, site_sc in enumerate(supercell): 579 d = supercell.lattice.get_distance_and_image(fcoords, site_sc.frac_coords)[0] 580 if d < 1e-5: 581 break 582 else: 583 raise RuntimeError(f"could not find a match for: {lines[iline]}") 584 585 for k in range(3): 586 fc[iat, isc, k] = [float(sp) for sp in lines[iline + 2 + k].split()] 587 iline += 5 588 return fc 589 590 591def parse_tdep_lotosplitting(filepath): 592 """ 593 Parses the lotosplitting file produced by TDEP and transforms them in 594 the phonopy format for Born effective charges and dielectric tensor. 595 596 Args: 597 filepath: path to the lotosplitting file. 598 599 Returns: 600 a tuple with dielectric tensor and Born effective charges. 601 """ 602 with open(filepath) as f: 603 lines = f.readlines() 604 605 values = [] 606 for l in lines: 607 if not l.strip(): 608 continue 609 values.append([float(v) for v in l.split()[:3]]) 610 611 if len(values) % 3 != 0: 612 raise RuntimeError(f"The data parsed has an unexpected shape: {np.shape(values)}") 613 614 eps = np.array(values[:3]) 615 natoms = len(values) // 3 - 1 616 born = np.reshape(values[3:], (natoms, 3, 3)) 617 618 return eps, born 619 620 621def write_tdep_lotosplitting(eps, born, filepath="infile.lotosplitting", fmt="%14.10f"): 622 """ 623 Writes an lotosplitting file starting from arrays containing dielectric tensor and 624 Born effective charges. 625 626 Args: 627 eps: a 3x3 array with the dielectric tensor. 628 born: an array with the Born effective charges. 629 filepath: the path where the lotosplitting file should be written. 630 fmt: the format for the float numbers. 631 """ 632 values = np.concatenate((eps[None, :, :], born), axis=0) 633 values = values.reshape((3 * len(values), 3)) 634 np.savetxt(filepath, values, fmt=fmt) 635 636 637def born_to_lotosplitting(born, lotosplitting_path="infile.lotosplitting"): 638 """ 639 Converted of a file from the BORN file produced from phonopy to the lotosplitting 640 file used by TDEP. 641 642 Args: 643 born: a dictionary with "dielectric" and "born" keywords as obtained from the nac_params 644 in phonopy. Notice that the "factor" attribute is not taken into account, so the 645 values should be in default phonopy units. 646 lotosplitting_path: the path where the lotosplitting file should be written. 647 """ 648 649 eps = born["dielectric"] 650 becs = born["born"] 651 write_tdep_lotosplitting(eps, becs, lotosplitting_path) 652 653 654@requires(Phonopy, "phonopy not installed!") 655def write_BORN(primitive, borns, epsilon, filename="BORN", symmetrize_tensors=False): 656 """ 657 Helper function imported from phonopy.file_IO. 658 Contrarily to the original, it does not symmetrize the tensor. 659 """ 660 lines = get_BORN_lines(primitive, borns, epsilon, symmetrize_tensors=symmetrize_tensors) 661 with open(filename, 'w') as w: 662 w.write('\n'.join(lines)) 663 664 665@requires(Phonopy, "phonopy not installed!") 666def get_BORN_lines(unitcell, borns, epsilon, 667 factor=None, 668 primitive_matrix=None, 669 supercell_matrix=None, 670 symprec=1e-5, symmetrize_tensors=False): 671 """ 672 Helper function imported from phonopy.file_IO that exposes the 673 option of not symmetrizing the tensor. 674 """ 675 from phonopy.structure.symmetry import elaborate_borns_and_epsilon 676 borns, epsilon, atom_indices = elaborate_borns_and_epsilon( 677 unitcell, borns, epsilon, symmetrize_tensors=symmetrize_tensors, 678 primitive_matrix=primitive_matrix, 679 supercell_matrix=supercell_matrix, 680 symprec=symprec) 681 682 text = "# epsilon and Z* of atoms " 683 text += ' '.join(["%d" % n for n in atom_indices + 1]) 684 lines = [text, ] 685 lines.append(("%13.8f " * 9) % tuple(epsilon.flatten())) 686 for z in borns: 687 lines.append(("%13.8f " * 9) % tuple(z.flatten())) 688 return lines 689