Merge pull request #140 from issp-center-dev/interaction_file

[model] interaction=file

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 <https://sites.google.com/view/kazuma7k6r>`_.
 **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
 ------------------------------------

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'))
+        print("\n    Written to {0}".format(p["model"]["seedname"]+'.h5'))

src/dcore/program_options.py

@@ -72,14 +72,15 @@ 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",
                       "U (Diagonal Coulomb pot.), U\' (Off-diagonal Coulomb pot.) and J (Hund coupling) (See below).")
     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)")