Merge pull request #11 from Green-Phys/ir_upgrade

switch mbtools to new IR grids while also keeping support for old format

green_mbtools/pesto/ir.py

@@ -4,14 +4,30 @@
 import h5py
 import os
 
-'''
-Fourier transform between imaginary time and Matsubara frequency using
-the intermediate representation (IR)
-'''
-
 
 class IR_factory(object):
-    def __init__(self, beta, ir_file=None):
+    """
+    Fourier transform between imaginary time and Matsubara frequency using
+    the intermediate representation (IR) grids.
+    """
+
+    def __init__(self, beta, ir_file=None, legacy_ir=False):
+        """Initialize IR basis class for post processing of Green data.
+
+        Parameters
+        ----------
+        beta : float
+            inverse temperature of the calculation
+        ir_file : string, required
+            IR-grid file, by default None
+        legacy : bool, optional
+            toggle legacy format for IR grid file, by default False
+
+        Raises
+        ------
+        ValueError
+            when a valid IR-grid file is not provided
+        """
 
         if ir_file is None:
             raise ValueError(
@@ -21,8 +37,12 @@ def __init__(self, beta, ir_file=None):
 
         self.beta = beta
         self.ir_file = ir_file
+        if legacy_ir:
+            self.read_IR_matrices = legacy_read_IR_matrices
+        else:
+            self.read_IR_matrices = new_read_IR_matrices
         self.tau_mesh, self.wsample, self.Ttc, self.Tcn, \
-            self.Tnc, self.Tct = read_IR_matrices(
+            self.Tnc, self.Tct = self.read_IR_matrices(
                 os.path.abspath(ir_file), self.beta
             )
         self.nts = self.tau_mesh.shape[0]
@@ -36,7 +56,7 @@ def update(self, beta=None, ir_file=None):
         if beta is not None:
             self.beta = beta
         self.tau_mesh, self.wsample, self.Ttc, self.Tcn, \
-            self.Tnc, self.Tct = read_IR_matrices(self.ir_file, self.beta)
+            self.Tnc, self.Tct = self.read_IR_matrices(self.ir_file, self.beta)
         self.nts = self.tau_mesh.shape[0]
         self.nw = self.wsample.shape[0]
 
@@ -76,7 +96,7 @@ def tauf_to_wb(self, X_t):
         E.g., P0(tau) -> P0(i Omega).
         """
 
-        _, wsample_bose, Ttc_b, _, Tnc_b, Tct_b = read_IR_matrices(
+        _, wsample_bose, Ttc_b, _, Tnc_b, Tct_b = self.read_IR_matrices(
                 os.path.abspath(self.ir_file), self.beta, ptype='bose'
             )
         fir = h5py.File(self.ir_file, 'r')
@@ -101,7 +121,7 @@ def wb_to_tauf(self, X_wb):
         """Transform quantity from bosonic frequency grid to fermionic
         tau-grid, e.g., P0(i Omega) -> P0(tau)
         """
-        _, wsample_bose, _, Tcn_b, _, _ = read_IR_matrices(
+        _, wsample_bose, _, Tcn_b, _, _ = self.read_IR_matrices(
                 os.path.abspath(self.ir_file), self.beta, ptype='bose'
             )
         nw_b = len(wsample_bose)
@@ -141,7 +161,36 @@ def tau_to_w_other(self, X_t, wsample):
         return X_w
 
 
-def read_IR_matrices(ir_path, beta, ptype='fermi'):
+def new_read_IR_matrices(ir_path, beta, ptype='fermi'):
+    ir = h5py.File(ir_path, 'r')
+    wsample = ir[ptype + "/ngrid"][()]
+    xsample = ir[ptype + "/xgrid"][()]
+
+    Ttc_minus1 = ir[ptype + "/u1l_neg"][()]
+    Ttc_tmp = ir[ptype + "/uxl"][()]
+    Ttc_1 = ir[ptype + "/u1l_pos"][()]
+    Ttc = np.zeros((Ttc_tmp.shape[0]+2, Ttc_tmp.shape[1]))
+    Ttc[0], Ttc[1:-1], Ttc[-1] = Ttc_minus1, Ttc_tmp, Ttc_1
+    Tnc = ir[ptype + "/uwl"][()]
+    ir.close()
+
+    if ptype == 'fermi':
+        zeta = 1
+    else:
+        zeta = 0
+    wsample = (2*wsample + zeta) * np.pi / beta
+    tau_mesh = np.zeros(xsample.shape[0]+2)
+    tau_mesh[0], tau_mesh[1:-1], tau_mesh[-1] = 0, (xsample+1)*beta/2.0, beta
+
+    Ttc *= np.sqrt(2.0/beta)
+    Tnc *= np.sqrt(beta)
+    Tct = np.linalg.inv(Ttc[1:-1])
+    Tcn = np.linalg.inv(Tnc)
+
+    return tau_mesh, wsample, Ttc, Tcn, Tnc, Tct
+
+
+def legacy_read_IR_matrices(ir_path, beta, ptype='fermi'):
     ir = h5py.File(ir_path, 'r')
     wsample = ir[ptype + "/wsample"][()]
     xsample = ir[ptype + "/xsample"][()]

green_mbtools/pesto/mb.py

@@ -14,7 +14,7 @@ class MB_post(object):
     """Many-body analysis class"""
     def __init__(
         self, fock, sigma=None, mu=None, gtau=None, S=None, kmesh=None,
-        beta=None, ir_file=None
+        beta=None, ir_file=None, legacy_ir=False
     ):
         """ Initialization """
         # Public instance variables
@@ -36,6 +36,7 @@ def __init__(
         self._ir_file = None
         self._beta = None
         self._mu = None
+        self.legacy_ir = legacy_ir
 
         """Setup"""
         if fock.ndim == 4:
@@ -114,7 +115,7 @@ def beta(self, value):
         print("Updated beta = {}".format(value))
         self._beta = value
         if self.ir is None:
-            self.ir = IR_factory(self.beta, self.ir_file)
+            self.ir = IR_factory(self.beta, self.ir_file, self.legacy_ir)
         else:
             self.ir.update(self.beta, self.ir_file)
 
@@ -132,7 +133,7 @@ def ir_file(self, value):
         print("Setting up IR grid for {}".format(value))
         self._ir_file = value
         if self.ir is None:
-            self.ir = IR_factory(self.beta, self.ir_file)
+            self.ir = IR_factory(self.beta, self.ir_file, self.legacy_ir)
         else:
             self.ir.update(self.beta, self.ir_file)
         self._nts = self.ir.nts
@@ -293,13 +294,15 @@ def AC_nevanlinna(
     ):
         """
         Analytical continuation using Nevanlinna interpolation
-        :param gtau_orth: imaginary time Green's function in the orthogonal basis, will be obtained from curren self.gtau if None
+        :param gtau_orth: imaginary time Green's function in the orthogonal
+            basis, will be obtained from curren self.gtau if None
         :param n_real: number of real frequency points
         :param w_min: smallest value on real frequency grid
         :param w_max: largest value on real frequency grid
         :param eta: broadening parameter
         :param outdir: [DEPRECATED]
-        :return: real frequency grid along with spectral function for a given Green's function
+        :return: real frequency grid along with spectral function for a given
+            Green's function
         """
         if gtau_orth is None:
             gtau_orth = orth.sao_orth(
@@ -375,7 +378,28 @@ def to_full_bz(X, conj_list, ir_list, bz_index, k_ind):
     return Y
 
 
-def initialize_MB_post(sim_path=None, input_path=None, ir_file=None):
+def initialize_MB_post(sim_path, input_path, ir_file, legacy_ir=False):
+    """Automatically reads the input and sim files, and initializes
+    a MB_post object for post processing of GREEN data.
+
+    Parameters
+    ----------
+    sim_path : string
+        results or sim file from the GREEN calculation
+    input_path : string
+        input file used in the GREEN calculation
+    ir_file : string
+        IR-grid file used in the GREEN calculation
+    legacy_ir : bool, optional
+        Toggle old format IR-grid file, by default False
+
+    Returns
+    -------
+    MB_post instance
+        object that contains information about the input and simulation,
+        and allows further post processing such as wannier interpolation,
+        analytic cotninuation, orthogonalization, etc.
+    """
     import h5py
     f = h5py.File(sim_path, 'r')
     it = f["iter"][()]
@@ -398,15 +422,20 @@ def initialize_MB_post(sim_path=None, input_path=None, ir_file=None):
     conj_list = f["grid/conj_list"][()]
     f.close()
 
-    """ All k-dependent matrices should lie on a full Monkhorst-Pack grid. """
+    """
+    All k-dependent matrices should lie on a full Monkhorst-Pack grid.
+    """
     F = to_full_bz(Fr, conj_list, ir_list, index, 1)
     S = to_full_bz(Sr, conj_list, ir_list, index, 1)
     Sigma = to_full_bz(Sigmar, conj_list, ir_list, index, 2)
     G = to_full_bz(Gr, conj_list, ir_list, index, 2)
 
-    """ Results from correlated methods """
+    """
+    Results from correlated methods
+    """
 
     # Standard way to initialize
     return MB_post(
-        fock=F, sigma=Sigma, mu=mu, gtau=G, S=S, beta=beta, ir_file=ir_file
+        fock=F, sigma=Sigma, mu=mu, gtau=G, S=S, beta=beta, ir_file=ir_file,
+        legacy_ir=legacy_ir
     )

tests/conftest.py

@@ -52,7 +52,8 @@ def mbo(data_path):
 
     # Standard way to initialize
     mbobj = mb.MB_post(
-        fock=F, sigma=Sigma, mu=mu, gtau=G, S=S, beta=1000, ir_file=irf
+        fock=F, sigma=Sigma, mu=mu, gtau=G, S=S, beta=1000,
+        ir_file=irf, legacy_ir=True
     )
 
     return mbobj

tests/ir_test.py

@@ -1,11 +1,14 @@
 from numpy.linalg import norm
+from green_mbtools.pesto.ir import IR_factory
+import h5py
+import pytest
 
 #
 # Test function for IR-Transform
 #
 
 
-def test_ir_transform(mbo):
+def test_ir_fourier_transform(mbo):
     """Transfrom from tau to iw and back.
     """
 
@@ -28,3 +31,26 @@ def test_ir_transform(mbo):
     # Transform back to G_tau
     gtau2 = my_ir.w_to_tau(G_iw)
     assert norm(gtau2 - gtau) < 1e-10
+
+
+def test_new_ir_fourier_transform(data_path):
+    """Tests 1) initialization of new-format IR grids, and
+    2) trnasforms G from tau to iw, and back to tau to check
+    functionality and accuracy.
+    """
+
+    # init IR handler
+    beta = 100
+    ir_file = pytest.test_data_dir + '/ir_grid/1e4_new.h5'
+    myir = IR_factory(beta, ir_file)
+
+    # load new-format GF2 data
+    gf2_path = pytest.test_data_dir + '/H2_GF2_new_ir/data.h5'
+    fdata = h5py.File(gf2_path, 'r')
+    g_tau = fdata['G_tau'][()]
+
+    # transform
+    g_iw = myir.tau_to_w(g_tau)
+    g_tau2 = myir.w_to_tau(g_iw)
+
+    assert norm(g_tau2 - g_tau) < 1e-10

tests/mbpost_test.py

@@ -16,7 +16,7 @@ def test_initMBpost(mbo, data_path):
     input_path = data_dir + '/H2_GW/input.h5'
     ir_path = data_dir + '/ir_grid/1e4_104.h5'
     mbo_new = mb.initialize_MB_post(
-        sim_path, input_path, ir_path
+        sim_path, input_path, ir_path, legacy_ir=True
     )
 
     # compare G