diff --git a/doc/reference/interaction.rst b/doc/reference/interaction.rst index 6d6ee6d5..f0e47b4a 100644 --- a/doc/reference/interaction.rst +++ b/doc/reference/interaction.rst @@ -14,6 +14,7 @@ Model Hamiltonian is defined as where .. math:: + :label: H_int {\hat H}_{\rm int} = \frac{1}{2} \sum_{i, \alpha \beta \gamma \delta,\sigma \sigma'} @@ -136,6 +137,58 @@ The basis can be specified by ``slater_basis`` parameter. See :doc:`slater_basis Use the output by `RESPACK `_. **Under construction.** +``interaction = file`` +------------------------------------ + +One can set the U tensor :math:`U^{i}_{\alpha \beta \gamma \delta}` from an external file. +The input parameters for this are + +:: + + interaction = file + interaction_file = ['file1.npy', 'file2.npy',] + +The filenames for all inequivalent shells are listed in ``interaction_file`` parameter. +The suppored formats are + +- **NumPy .npy format (binary file)**: Four dimensional NumPy array corresponding to :math:`U^{i}_{\alpha \beta \gamma \delta}` for shell :math:`i` is save by ``numpy.save`` command. + +- **text file**: :math:`U^{i}_{\alpha \beta \gamma \delta}` for each shell :math:`i` is written in a one dimansional manner as follows (example for ``norb=3`` without spin-orbit coupling) + + :: + + Re(U[0,0,0,0]) + Re(U[0,0,0,1]) + Re(U[0,0,0,2]) + Re(U[0,0,1,0]) + ... + + or + + :: + + Re(U[0,0,0,0]) Im(U[0,0,0,0]) + Re(U[0,0,0,1]) Im(U[0,0,0,1]) + Re(U[0,0,0,2]) Im(U[0,0,0,2]) + Re(U[0,0,1,0]) Im(U[0,0,1,0]) + ... + + No blank line should be put between consecutive elements. Lines beginning with `#` are skipped. + +The definition of the U tensor depends on the ``spin_orbit`` parameter. +For ``spin_orbit=False``, :math:`U^{i}_{\alpha \beta \gamma \delta}` is defined by Eq. :eq:`H_int`. The shape of the U tensor is ``(norb, norb, norb, norb)``. + +For ``spin_orbit=True``, on the other hand, :math:`U^{i}_{\alpha \beta \gamma \delta}` is defined by + +.. math:: + {\hat H}_{\rm int} = \frac{1}{2} + \sum_{i, \alpha \beta \gamma \delta} + U^{i}_{\alpha \beta \gamma \delta} + c_{i \alpha}^\dagger c_{i \beta}^\dagger c_{i \delta} c_{i \gamma}. + +In this case, the indices :math:`\alpha, \beta, \gamma, \delta` include spin. +The shape of the U tensor is therefore ``(2*norb, 2*norb, 2*norb, 2*norb)``. + ``density_density`` option ------------------------------------ diff --git a/src/dcore/interaction.py b/src/dcore/interaction.py index d436e7f7..c340c822 100644 --- a/src/dcore/interaction.py +++ b/src/dcore/interaction.py @@ -290,7 +290,7 @@ def _generate_umat_slater(p: Dict, l_sh: List[int], f_sh: List[numpy.ndarray]): def _generate_umat_slater_uj(p: Dict): assert isinstance(p, Dict) - _check_parameters(p['model'], required=['slater_uj'], unused=['slater_f', 'kanamori']) + _check_parameters(p['model'], required=['slater_uj'], unused=['slater_f', 'kanamori', 'interaction_file']) def _U_J_to_F(_l, _u, _j): if _l == 0: @@ -310,7 +310,7 @@ def _U_J_to_F(_l, _u, _j): def _generate_umat_slater_f(p: Dict): assert isinstance(p, Dict) - _check_parameters(p['model'], required=['slater_f'], unused=['slater_uj', 'kanamori']) + _check_parameters(p['model'], required=['slater_f'], unused=['slater_uj', 'kanamori', 'interaction_file']) # parse slater parameters nsh = p['model']['n_inequiv_shells'] @@ -329,7 +329,7 @@ def _generate_umat_slater_f(p: Dict): def _generate_umat_respack(p: Dict): assert isinstance(p, Dict) - _check_parameters(p['model'], required=[], unused=['kanamori', 'slater_f', 'slater_uj']) + _check_parameters(p['model'], required=[], unused=['kanamori', 'slater_f', 'slater_uj', 'interaction_file']) nsh = p['model']['n_inequiv_shells'] norb = p['model']['norb_inequiv_sh'] @@ -382,6 +382,68 @@ def _generate_umat_respack(p: Dict): return u_mat_so_sh +def _parse_filenames(input_str: str, name, nsh: int): + try: + files = ast.literal_eval(input_str) + assert isinstance(files, list), f"interpreted as {type(files)}. List is required!" + assert len(files) == nsh, f"length={len(files)} must be {nsh} (# of inequivalent shells)" + except Exception as e: + print(f"\nError: {name} = {repr(input_str)}", file=sys.stderr) + print(e, file=sys.stderr) + exit(1) + print(f" {name} = {repr(files)}") + return files + + +def _generate_umat_file(p: Dict): + assert isinstance(p, Dict) + _check_parameters(p['model'], required=['interaction_file'], unused=['slater_f', 'slater_uj', 'kanamori']) + + nsh = p["model"]['n_inequiv_shells'] + norb_sh = p['model']['norb_inequiv_sh'] + spin_orbit = p["model"]["spin_orbit"] + + # str -> List[str] + filenames = _parse_filenames(p["model"]["interaction_file"], name='interaction_file', nsh=nsh) + + u_mat_so_sh = [] + print(f"\n Reading U-matrix from file(s)") + for ish, file in enumerate(filenames): + norb = norb_sh[ish] + print(f" ish={ish} : '{file}'") + + u_shape = (2*norb, 2*norb, 2*norb, 2*norb) if spin_orbit else (norb, norb, norb, norb) + + # read file + try: + if os.path.splitext(file)[1] == ".npy": + umat = numpy.load(file) + assert umat.shape == u_shape, f"inconsistent shape: require {u_shape}, but {umat.shape} is given." + else: + umat_1d = numpy.loadtxt(file) + if umat_1d.shape == (numpy.prod(u_shape),): # real U + umat = umat_1d.reshape(u_shape) + elif umat_1d.shape == (numpy.prod(u_shape), 2): # complex U + # float (norb^4, 2) -> complex (norb^4, 1) -> complex (norb,norb,norb,norb) + umat = umat_1d.view(complex).reshape(u_shape) + else: + raise Exception(f"inconsistent shape: require {(numpy.prod(u_shape),)} for real U or {(numpy.prod(u_shape), 2)} for complex U, but {umat_1d.shape} is given.") + except Exception as e: + print(f"\nError in reading file '{file}' for ish={ish}", file=sys.stderr) + print(e, file=sys.stderr) + exit(1) + + if spin_orbit: + umat_so = umat + else: + umat_so = to_spin_full_U_matrix(umat) + assert umat_so.shape == (2*norb, 2*norb, 2*norb, 2*norb) + + u_mat_so_sh.append(umat_so) + + return u_mat_so_sh + + def generate_umat(p: Dict): """ Add U-matrix block (Tentative) @@ -401,6 +463,7 @@ def generate_umat(p: Dict): 'slater_uj': _generate_umat_slater_uj, 'slater_f': _generate_umat_slater_f, 'respack': _generate_umat_respack, + 'file': _generate_umat_file, }.get(interaction) if func_umat is None: sys.exit(f"Error ! Invalid interaction : {interaction}") @@ -429,4 +492,4 @@ def generate_umat(p: Dict): f.create_group("DCore") f["DCore"]["Umat"] = u_mat_so_sh - print("\n Written to {0}".format(p["model"]["seedname"]+'.h5')) \ No newline at end of file + print("\n Written to {0}".format(p["model"]["seedname"]+'.h5')) diff --git a/src/dcore/program_options.py b/src/dcore/program_options.py index f6fba8b4..c13e7ecd 100644 --- a/src/dcore/program_options.py +++ b/src/dcore/program_options.py @@ -72,7 +72,7 @@ def create_parser(target_sections=None): parser.add_option("model", "nk2", int, 0, "Number of *k* along b_2 (for lattice = wannier90, external)") parser.add_option("model", "spin_orbit", bool, False, "Whether the spin-orbit case.") parser.add_option("model", "interaction", str, "kanamori", - 'Chosen from "slater_uj", "slater_f", "kanamori", "respack" (See below)') + 'Chosen from "slater_uj", "slater_f", "kanamori", "respack", "file" (See below)') parser.add_option("model", "density_density", bool, False, "If true, only the density-density part of the interaction is used (See below).") parser.add_option("model", "kanamori", str, "None", @@ -80,6 +80,7 @@ def create_parser(target_sections=None): parser.add_option("model", "slater_f", str, "None", "Angular momentum, Slater integrals F (See below).") parser.add_option("model", "slater_uj", str, "None", "Angular momentum, Slater integrals in U and J (See below).") parser.add_option("model", "slater_basis", str, "cubic", "Basis of the Slater interaction (See below).") + parser.add_option("model", "interaction_file", str, "None", "list of filenames that contain U tensor for each inequivalent shell") parser.add_option("model", "local_potential_matrix", str, "None", "dict of {ish: 'filename'} to specify local potential matrix of ish-th shell") parser.add_option("model", "local_potential_factor", str, "1.0", "Prefactors to the local potential matrix (float or list with len=ncor)")