More grids for Dynamic representations

dyson/__init__.py

@@ -71,6 +71,8 @@
        self-energy.
    * - :class:`~dyson.solvers.dynamic.cpgf.CPGF`
      - Chebyshev polynomial moments of the dynamic Green's function.
+   * - :class:`~dyson.solvers.dynamic.direct.Direct`
+     - Dynamic self-energy and initial Green's function, in frequency space.
 
 For a full accounting of the inputs and their types, please see the documentation for each solver.
 
@@ -135,5 +137,6 @@
     DensityRelaxation,
     CorrectionVector,
     CPGF,
+    Direct,
 )
 from dyson.expressions import HF, CCSD, FCI, ADC2, ADC2x, Hamiltonian

dyson/grids/__init__.py

@@ -11,7 +11,73 @@
 
     grid
     frequency
+    time
+    fourier
+    pade
+    util
 """
 
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, TypeVar
+
+from dyson import numpy as np
+from dyson.grids.grid import BaseGrid
 from dyson.grids.frequency import RealFrequencyGrid, GridRF
 from dyson.grids.frequency import ImaginaryFrequencyGrid, GridIF
+from dyson.grids.time import RealTimeGrid, GridRT
+from dyson.grids.time import ImaginaryTimeGrid, GridIT
+from dyson.grids.fourier import fourier_transform_imag
+from dyson.grids.pade import analytic_continuation_freq_pade
+from dyson.grids.util import are_dual, are_equal
+
+if TYPE_CHECKING:
+    from typing import Any
+
+    from dyson.representations import Dynamic
+
+GridSrcT = TypeVar("GridSrcT", bound=BaseGrid)
+GridDstT = TypeVar("GridDstT", bound=BaseGrid)
+
+
+def transform(dynamic: Dynamic[GridSrcT], grid: GridDstT, **kwargs: Any) -> Dynamic[GridDstT]:
+    """Transform a dynamic quantity to a new grid using either FFT or AC.
+
+    Currently available transformations are:
+
+    .. code-block:: bash
+
+                   AC
+               ─────────>
+        GridRF <───────── GridIF
+                   AC
+                           │ ^
+                           │ │
+                      IFFT │ │ FFT
+                           │ │
+                           v │
+
+                         GridIT
+
+    Args:
+        dynamic: The dynamic quantity to transform.
+        grid: The grid to transform to.
+        **kwargs: Additional keyword arguments passed to the transformation function.
+
+    Returns:
+        The transformed dynamic quantity.
+
+    Raises:
+        NotImplementedError: If the transformation is not implemented.
+    """
+    if isinstance(dynamic.grid, GridIT) and isinstance(grid, GridIF):
+        return fourier_transform_imag(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridIF) and isinstance(grid, GridIT):
+        return fourier_transform_imag(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridIF) and isinstance(grid, GridRF):
+        return analytic_continuation_freq_pade(dynamic, grid, **kwargs)  # type: ignore
+    if isinstance(dynamic.grid, GridRF) and isinstance(grid, GridIF):
+        return analytic_continuation_freq_pade(dynamic, grid, **kwargs)  # type: ignore
+    raise NotImplementedError(
+        f"transformation between {type(dynamic.grid)} and {type(grid)} not implemented"
+    )

dyson/grids/fourier.py

@@ -0,0 +1,81 @@
+"""Fourier transformation."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from dyson import numpy as np
+from dyson import util
+from dyson.grids.util import are_dual
+from dyson.representations.enums import Component, Reduction
+
+if TYPE_CHECKING:
+    from dyson.grids.grid import BaseImaginaryGrid
+    from dyson.representations.dynamic import Dynamic
+
+
+def fourier_transform_imag(
+    greens_function: Dynamic[BaseImaginaryGrid], grid: BaseImaginaryGrid
+) -> Dynamic[BaseImaginaryGrid]:
+    """Fourier transform between imaginary frequency and imaginary time grids.
+
+    Args:
+        greens_function: Dynamic quantity in imaginary frequency or time domain.
+        grid: Target grid for the Fourier transform.
+
+    Returns:
+        Dynamic quantity in the target domain.
+    """
+    grid_in, grid_out = greens_function.grid, grid
+    if not are_dual(grid_in, grid_out):
+        raise ValueError("the two grids must be dual to each other.")
+    if grid_in.reality or grid_out.reality:
+        raise ValueError("only imaginary frequency and imaginary time grids is supported.")
+    if greens_function.component != Component.FULL:
+        raise ValueError("only full component is supported.")
+    if greens_function.reduction == Reduction.TRACE:
+        raise ValueError("traced reduction is not supported.")
+    forward = grid_in.domain == "time"
+
+    # Setup based on direction
+    beta = grid_in.beta
+    if forward:
+        sign = 1
+        norm = beta
+        fft = np.fft.ifft
+    else:
+        sign = -1
+        norm = 2 / beta
+        fft = np.fft.fft
+
+    # Get the shifts
+    shifts = (
+        np.exp(1.0j * sign * np.min(grid_out.points) * grid_in.points),
+        np.exp(1.0j * sign * np.min(grid_in.points) * (grid_out.points - np.min(grid_out.points))),
+    )
+
+    # Get the input array
+    array = greens_function.array.copy()
+
+    # Perform the Fourier transform
+    array = util.einsum("w...,w->w...", array, shifts[0])
+    array = fft(array, max(len(grid_in), len(grid_out)), axis=0)[: len(grid_out)]
+    array = util.einsum("w...,w->w...", array, shifts[1])
+
+    # Normalise
+    array *= norm
+
+    # Hermitise for imaginary frequency to imaginary time transform
+    if not forward:
+        if greens_function.reduction == Reduction.NONE:
+            array = 0.5 * (array + array.swapaxes(1, 2).conj())
+        else:
+            array = 0.5 * (array + array.conj())
+
+    return greens_function.__class__(
+        grid_out,
+        array,
+        reduction=greens_function.reduction,
+        ordering=greens_function.ordering,
+        hermitian=greens_function.hermitian,
+    )