Repository: ochsenfeld-lab/adaptive_sampling (language: python)

Issue:
Switch from multiply to matmul for speedup
Rewrote multiplications and summations in the mbar loop using the optimized matmul for a significant speedup even on CPUs

Also, float64 seems to be necessary for convergence as float32 can diverge for huge amounts of data.


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 672ca43..71bf8f2 100644
--- a/README.md
+++ b/README.md
@@ -84,7 +84,7 @@ the_bias = eABF(
 for md_step in range(steps):
     # propagate langevin dynamics and calc forces 
     ... 
-    bias_force = eABF.step_bias(write_output=True, write_traj=True)
+    bias_force = the_bias.step_bias(write_output=True, write_traj=True)
     the_md.forces += bias_force
     ...
     # finish md_step
@@ -123,6 +123,7 @@ weights = mbar.run_mbar(
     conv=1.0e-7,
     conv_errvec=1.0,
     outfreq=100,
+    device='cpu',
 )
 
 pmf, rho = mbar.pmf_from_weights(grid, cv[indices], weights, equil_temp=300.0)
diff --git a/adaptive_sampling/colvars/__init__.py b/adaptive_sampling/colvars/__init__.py
index 12449a4..d1dd66d 100644
--- a/adaptive_sampling/colvars/__init__.py
+++ b/adaptive_sampling/colvars/__init__.py
@@ -1,18 +1,7 @@
-from .colvars import CV
-from .kearsley import Kearsley
+from .colvars import * 
 from .utils import *
 
 __all__ = [
     "colvars",
-    "kearsley",
     "utils",
 ]
-from .colvars import *
-from .kearsley import * 
-from .utils import *
-
-__all__ = [
-    "colvars",
-    "kearsley",
-    "utils",
-]
\ No newline at end of file
diff --git a/adaptive_sampling/colvars/colvars.py b/adaptive_sampling/colvars/colvars.py
index c1685e8..f4dcf99 100644
--- a/adaptive_sampling/colvars/colvars.py
+++ b/adaptive_sampling/colvars/colvars.py
@@ -6,7 +6,6 @@
 from ..interface.sampling_data import MDInterface
 from ..units import *
 from .utils import *
-from .kearsley import Kearsley
 
 
 class CV:
@@ -336,8 +335,6 @@ def rmsd(
 
         if method.lower() == "kabsch":
             self.cv = kabsch_rmsd(self.coords, self.reference, indices=atom_indices)
-        elif method.lower() == "kearsley":
-            self.cv = Kearsley().fit(self.coords, self.reference, indices=atom_indices)
         else:  # 'quaternion':
             self.cv = quaternion_rmsd(self.coords, self.reference, indices=atom_indices)
 
@@ -349,288 +346,6 @@ def rmsd(
 
         return float(self.cv)
 
-    def path_progression(
-        self,
-        cv_def: Union[str, list],
-        n_interpol: int = 2,
-        method: str = "internal",
-    ) -> float:
-        """progression along path
-
-        see: Branduardui et al., J. Chem. Phys. (2007); https://doi.org/10.1063/1.2432340
-             Leines et. al., Phys. Rev. Lett. (2012), https://link.aps.org/doi/10.1103/PhysRevLett.109.020601
-
-        TODO: Adaptive path updates to converge to MFEP 
-
-        Args:
-            cv_def: path to xyz file with reference structure
-                    definition: 'path to xyz' or
-                                ['path to xyz', [atom indices]]
-            n_interpol: number of interpolated images to between two nodes of the reference path
-            method: 'gpath' for geometrical path definition based on selected internal distances (Leines et. al.),
-                    'internal' for using RMSD of selected internal distances, 
-                    'quaternion', 'kabsch' or 'kearsley' for RMSD based on optimal alignment of cartesian coordinates,
-
-        Returns:
-            cv: path collective variable
-        """
-        self.update_coords()
-
-        if isinstance(cv_def, list):
-            atom_indices = cv_def[1]
-            cv_def = cv_def[0]
-        else:
-            atom_indices = None
-
-        if self.reference == None:
-            images = read_traj(cv_def)
-            self.reference = interpolate_coordinates(images, n_interpol=n_interpol)
-            if len(self.reference[0]) != len(self.coords):
-                raise ValueError(
-                    " >>> Error: number of cartesian coordinates has to match reference coordinates"
-                )
-
-            if method.lower() == 'gpath':
-                #self.close_index = find_closest_points(
-                #    self.coords, indices=atom_indices
-                #)
-                print(
-                    f"\n >>> Colvars Info: Geometrical Path CV defined by {len(self.reference)} nodes."
-                )
-            else:
-
-                print(
-                    f"\n >>> Colvars Info: Path CV defined by {len(self.reference)} nodes."
-                )
-                self.la = 0
-                if method.lower() == "internal":
-                    print(
-                        f" >>> Colvars Info: Using internal coordinates for calculation of RMSD."
-                    )
-                    if self.reference_internal == None:
-                        self.reference_internal = []
-                        for image in self.reference:
-                            self.reference_internal.append(
-                                get_internal_coords(
-                                    image, atom_indices,
-                                )
-                            )
-                    for i in range(1, len(self.reference_internal)):
-                        self.la += torch.linalg.norm(
-                            self.reference_internal[i] - self.reference_internal[i - 1]
-                        )
-
-                else:    
-                    print(
-                        f" >>> Colvars Info: Using {method} algorithm for calculation of RMSD."
-                    )
-                    for i in range(1, len(self.reference)):
-                        self.la += torch.linalg.norm(
-                            self.reference[i] - self.reference[i - 1]
-                        )
-
-                self.la /= float(len(self.reference))
-                self.la = 1.0 / torch.pow(self.la, 2)
-                print(
-                    f" >>> Colvars Info: Setting lambda parameter to {self.la} 1/Bohr^2."
-                )
-
-        M = len(self.reference)
-
-        if method.lower() == 'gpath':
-            # use geometrical path definition in internal coordinates
-            self.coords_internal = get_internal_coords(
-                self.coords, atom_indices,
-            )
-            if self.reference_internal == None:
-                self.reference_internal = []
-                for image in self.reference:
-                    self.reference_internal.append(
-                        get_internal_coords(
-                            image, atom_indices,
-                        )
-                    )
-
-            # find index of nearest image on path
-            rmsds = []
-            for i, image in enumerate(self.reference_internal):
-                rmsds.append(torch.linalg.norm(image - self.coords_internal))
-            rmsds = torch.stack(rmsds)
-            
-            nearest_image = torch.argmin(rmsds)
-            if int(nearest_image) == 0:
-                nearest_image += 1
-            elif int(nearest_image) == len(rmsds):
-                nearest_image -= 1
-
-            # compute path cv
-            v1 = self.reference_internal[nearest_image] - self.coords_internal
-            v2 = self.coords_internal - self.reference_internal[nearest_image - 1]
-            v3 = (
-                self.reference_internal[nearest_image + 1]
-                - self.reference_internal[nearest_image]
-            )
-
-            v1_n2 = torch.pow(torch.linalg.norm(v1), 2)
-            v2_n2 = torch.pow(torch.linalg.norm(v2), 2)
-            v3_n2 = torch.pow(torch.linalg.norm(v3), 2)
-
-            dot13 = torch.dot(v1, v3)
-            denom = 2.0 * M * v3_n2
-
-            self.cv = (
-                torch.sqrt(torch.pow(dot13, 2) - v3_n2 * (v1_n2 - v2_n2)) / denom
-            ) - ((dot13 - v3_n2) / denom)
-
-            if rmsds[nearest_image - 1] < rmsds[nearest_image + 1]:
-                self.cv = float(nearest_image) / float(M) - self.cv
-            else:
-                self.cv = float(nearest_image) / float(M) + self.cv
-
-        else:
-
-            if method.lower() == "internal":
-                self.coords_internal = get_internal_coords(
-                    self.coords, atom_indices,
-                )
-            
-            num = 0
-            denom = 0
-
-            for i, image in enumerate(self.reference):
-
-                if method.lower() == "internal":
-                    d = torch.linalg.norm(self.reference_internal[i] - self.coords_internal)
-                elif method.lower() == "kabsch":
-                    d = kabsch_rmsd(image, self.coords, indices=atom_indices)
-                elif method.lower() == "kearsley":
-                    d = Kearsley().fit(image, self.coords, indices=atom_indices)
-                elif method.lower() == 'quaternion':
-                    d = quaternion_rmsd(image, self.coords, indices=atom_indices)
-                else:
-                    raise ValueError(" >>> Error: invalid method for calculation of RMSD.")
-
-                exp = torch.exp(-self.la * d)
-                num += float(i) * exp
-                denom += exp
-
-            self.cv = num / denom
-            self.cv = self.cv / float(M)  # normalize to range(0,1)
-
-        # get forces
-        if self.requires_grad:
-            self.gradient = torch.autograd.grad(
-                self.cv, self.coords, allow_unused=True
-            )[0]
-            self.gradient = self.gradient.detach().numpy()
-
-        return float(self.cv)
-
-    def path_distance(
-        self,
-        cv_def: Union[str, list],
-        n_interpol: int = 2,
-        method: str = "quaternion",
-    ) -> float:
-        """distance from path 
-
-        see: Branduardui et al., J. Chem. Phys. (2007); https://doi.org/10.1063/1.2432340
-             
-        Args:
-            cv_def: path to xyz file with reference path
-                    definition: 'path to xyz' or
-                                ['path to xyz', [atom indices]]
-            n_interpol: number of interpolated images to between to nodes of path
-            method: 'internal' for using RMSD of selected internal distances, 
-                    'quaternion', 'kabsch' or 'kearsley' for using optimal alignment of cartesian coordinates,
-
-        Returns:
-            cv: distance from path
-        """
-        self.update_coords()
-
-        if isinstance(cv_def, list):
-            atom_indices = cv_def[1]
-            cv_def = cv_def[0]
-        else:
-            atom_indices = None
-
-        if self.reference == None:
-            images = read_traj(cv_def)
-            self.reference = interpolate_coordinates(images, n_interpol=n_interpol)
-            if len(self.reference[0]) != len(self.coords):
-                raise ValueError(
-                    " >>> Error: number of cartesian coordinates has to match reference coordinates"
-                )
-
-            print(
-                f"\n >>> Colvars Info: Distance to Path defined by {len(self.reference)} nodes."
-            )
-            self.la = 0
-            if method.lower() == "internal":
-                print(
-                    f" >>> Colvars Info: Using internal coordinates for calculation of RMSD."
-                )
-                if self.reference_internal == None:
-                    self.reference_internal = []
-                    for image in self.reference:
-                        self.reference_internal.append(
-                            get_internal_coords(
-                                image, atom_indices,
-                            )
-                        )
-                for i in range(1, len(self.reference_internal)):
-                    self.la += torch.linalg.norm(
-                        self.reference_internal[i] - self.reference_internal[i - 1]
-                    )
-
-            else:    
-                print(
-                    f" >>> Colvars Info: Using {method} algorithm for calculation of RMSD."
-                )
-                for i in range(1, len(self.reference)):
-                    self.la += torch.linalg.norm(
-                        self.reference[i] - self.reference[i - 1]
-                    )
-
-            self.la /= float(len(self.reference))
-            self.la = 1.0 / torch.pow(self.la, 2)
-            print(
-                f" >>> Colvars Info: Setting lambda parameter to {self.la} 1/Bohr^2."
-            )
-
-        if method.lower() == "internal":
-            self.coords_internal = get_internal_coords(
-                self.coords, atom_indices,
-            )
-
-        sum = 0
-        for i, image in enumerate(self.reference):
-
-            if method.lower() == "internal":
-                d = torch.linalg.norm(self.reference_internal[i] - self.coords_internal)
-            elif method.lower() == "kabsch":
-                d = kabsch_rmsd(image, self.coords, indices=atom_indices)
-            elif method.lower() == "kearsley":
-                d = Kearsley().fit(image, self.coords, indices=atom_indices)
-            elif method.lower() == 'quaternion':
-                d = quaternion_rmsd(image, self.coords, indices=atom_indices)
-            else:
-                raise ValueError(" >>> Error: invalid method for calculation of RMSD.")
-
-            sum += torch.exp(-self.la * d)
-
-        self.cv = - torch.log(sum) / self.la
-
-        # get forces
-        if self.requires_grad:
-            self.gradient = torch.autograd.grad(
-                self.cv, self.coords, allow_unused=True
-            )[0]
-            self.gradient = self.gradient.detach().numpy()
-
-        return float(self.cv)
-
     def get_cv(self, cv: str, atoms: list, **kwargs) -> Tuple[float, np.ndarray]:
         """get state of collective variable from cv definition of sampling_tools
 
@@ -672,18 +387,6 @@ def get_cv(self, cv: str, atoms: list, **kwargs) -> Tuple[float, np.ndarray]:
         elif cv.lower() == "rmsd":
             xi = self.rmsd(atoms)
             self.type = "distance"
-        elif cv.lower() == "path":
-            xi = self.path_progression(atoms, method='quaternion', **kwargs)
-            self.type = None
-        elif cv.lower() == "cv_path":
-            xi = self.path_progression(atoms, method='internal', **kwargs)
-            self.type = None
-        elif cv.lower() == "gpath":
-            xi = self.path_progression(atoms, method='gpath', **kwargs)
-            self.type = None
-        elif cv.lower() == "path_distance":
-            xi = self.path_distance(atoms, **kwargs)
-            self.type = None
         else:
             print(" >>> Error in CV: Unknown coordinate")
             sys.exit(1)
diff --git a/adaptive_sampling/colvars/kearsley.py b/adaptive_sampling/colvars/kearsley.py
deleted file mode 100644
index a4f1c0c..0000000
--- a/adaptive_sampling/colvars/kearsley.py
+++ /dev/null
@@ -1,202 +0,0 @@
-import torch
-from scipy.spatial.transform import Rotation
-
-
-class Kearsley:
-    """
-    A Kearsley transformation.
-    Algorithm to minimize the sum of the squared distances between two sets of points.
-    Author: Simon K. Kearsley
-    Paper: "On the orthogonal transformation used for structural comparisons"
-    https://doi.org/10.1107/S0108767388010128
-
-    Args:
-        rot: Rotation
-             A scipy.spatial.transform Rotation in 3 dimensions.
-        trans: ndarray, shape (3,)
-             The 3D translation.
-
-    Examples:
-
-    Given two sets of 3D points u and v:
-    >>> u
-    array([[0, 0, 0],
-           [0, 0, 1],
-           [0, 0, 2],
-           ...,
-           [9, 9, 7],
-           [9, 9, 8],
-           [9, 9, 9]])
-    >>> v
-    array([[ 30.50347534, -20.16089091,  -7.42752623],
-           [ 30.77704903, -21.02339348,  -7.27823201],
-           [ 31.3215374 , -21.99452332,  -7.15703548],
-           ...,
-           [ 42.05988643, -23.50924264, -15.59516355],
-           [ 42.27217891, -24.36478643, -15.59064995],
-           [ 42.66080502, -25.27318759, -15.386241  ]])
-    It is possible to calculate the rotation and translation that minimize the root mean squared deviation:
-    >>> from kearsley import Kearsley
-    >>> k = Kearsley()
-    >>> rmsd = k.fit(u, v)
-    >>> rmsd
-    0.10003430497284149
-    The rotation and translation are the attributes of the class:
-    >>> k.rot.as_matrix()
-    array([[ 0.05552838, -0.04405506, -0.99748471],
-           [ 0.91956342,  0.39147652,  0.03390061],
-           [ 0.38899835, -0.9191329 ,  0.06224948]])
-    >>> k.trans
-    array([ 30.46560753, -20.15086287,  -7.34422276])
-
-    Once fitted you can apply the transformation:
-    >>> v_transform = k.transform(v)
-    >>> v_transform
-    array([[ 0.08563846,  0.02807207,  0.01876202],
-           [-0.01009153, -0.0529479 ,  0.92722971],
-           [-0.05796549,  0.07167779,  2.03917659],
-           ...,
-           [ 9.0219524 ,  9.067236  ,  7.08333594],
-           [ 9.06692944,  8.9276801 ,  7.95255679],
-           [ 8.92463409,  8.93635832,  8.95139744]])
-    It is also possible to fit and transform with one command:
-    >>> v_transform, rmsd = k.fit_transform(u, v)
-    >>> rmsd
-    0.10003430497284149
-    >>> v_transform
-    array([[ 0.08563846,  0.02807207,  0.01876202],
-           [-0.01009153, -0.0529479 ,  0.92722971],
-           [-0.05796549,  0.07167779,  2.03917659],
-           ...,
-           [ 9.0219524 ,  9.067236  ,  7.08333594],
-           [ 9.06692944,  8.9276801 ,  7.95255679],
-           [ 8.92463409,  8.93635832,  8.95139744]])
-    """
-
-    def __init__(self, requires_grad: bool = True):
-        """The attributes are filled with a rotation of zero degrees and a translation of zero."""
-        self.requires_grad = requires_grad
-        self.rot = Rotation.from_quat([0, 0, 0, 1])
-        self.trans = torch.zeros(3)
-
-    def _kearsley_matrix(self, x: torch.tensor, y: torch.tensor) -> torch.tensor:
-        """
-        Calculates the Kearsley matrix.
-
-        Args:
-            x: shape (N, 3), Input array with 3D points, centered at zero.
-            y: shape (N, 3), Input array with 3D points, centered at zero.
-
-        Returns:
-            K: shape (4, 4), The Kearsley matrix.
-        """
-        # diff and sum quantities
-        d, s = x - y, x + y
-
-        # extract columns to simplify notation
-        d0, d1, d2 = d[:, 0], d[:, 1], d[:, 2]
-        s0, s1, s2 = s[:, 0], s[:, 1], s[:, 2]
-
-        # fill kearsley matrix
-        K = torch.zeros((4, 4))
-        K[0, 0] = torch.matmul(d0, d0) + torch.matmul(d1, d1) + torch.matmul(d2, d2)
-        K[1, 0] = torch.matmul(s1, d2) - torch.matmul(d1, s2)
-        K[2, 0] = torch.matmul(d0, s2) - torch.matmul(s0, d2)
-        K[3, 0] = torch.matmul(s0, d1) - torch.matmul(d0, s1)
-        K[1, 1] = torch.matmul(s1, s1) + torch.matmul(s2, s2) + torch.matmul(d0, d0)
-        K[2, 1] = torch.matmul(d0, d1) - torch.matmul(s0, s1)
-        K[3, 1] = torch.matmul(d0, d2) - torch.matmul(s0, s2)
-        K[2, 2] = torch.matmul(s0, s0) + torch.matmul(s2, s2) + torch.matmul(d1, d1)
-        K[3, 2] = torch.matmul(d1, d2) - torch.matmul(s1, s2)
-        K[3, 3] = torch.matmul(s0, s0) + torch.matmul(s1, s1) + torch.matmul(d2, d2)
-
-        return K
-
-    def fit(
-        self, u: torch.tensor, v: torch.tensor, indices: list = None
-    ) -> torch.tensor:
-        """
-        Calculates the rotation and translation that best fits both sets of points.
-
-        Args:
-            u: shape (3*N,), Input array with 3D points.
-            v: shape (3*N,), Input array with 3D points.
-            indices: list of indices that are included in calculation of rmsd
-
-        Returns:
-            rmsd: The root mean squared deviation.
-        """
-        u_r = torch.reshape(u, (int(len(u) / 3), 3))
-        v_r = torch.reshape(v, (int(len(v) / 3), 3))
-
-        if indices != None:
-            u_r = u_r[indices]
-            v_r = v_r[indices]
-
-        # centroids
-        centroid_u = u_r.mean(axis=0)
-        centroid_v = v_r.mean(axis=0)
-
-        # center both sets of points
-        x, y = u_r - centroid_u, v_r - centroid_v
-
-        # calculate Kearsley matrix
-        K = self._kearsley_matrix(x, y)
-
-        # diagonalize K
-        eig_vals, eig_vecs = torch.linalg.eigh(K)
-
-        # calculate rotation and translation
-        if self.requires_grad == False:
-            q = eig_vecs[:, 0]
-            q = torch.roll(q, 3)
-
-            # this would break autograd!
-            self.rot = Rotation.from_quat(q).inv()
-            self.trans = centroid_v - self.rot.inv().apply(centroid_u)
-
-        # calculate rmsd
-        eig_val = torch.abs(eig_vals[0])
-        rmsd = torch.sqrt(eig_val / u_r.size(dim=0))
-
-        return rmsd
-
-    def transform(self, u: torch.matmul) -> torch.matmul:
-        """
-        Transforms a list of 3D points with a rotation and a translation.
-
-        Args:
-            u: shape (N, 3), Input array with 3D points.
-
-        Returns:
-            t: shape (N, 3), Input points transformed.
-
-        Raises:
-            ValueError: If the input points don't have the correct shape.
-        """
-        if len(u.size(dim=0)) != 2:
-            raise ValueError("Input array must have 2 dimensions")
-
-        if u.sizre(dim=1) != 3:
-            raise ValueError("Input array must have 3 columns")
-
-        return self.rot.apply(u - self.trans)
-
-    def fit_transform(
-        self, u: torch.tensor, v: torch.tensor, indices: list = None
-    ) -> tuple:
-        """
-        Calculates the rotation and translation that best fits both sets of points and
-        applies the transformation to the second set.
-
-        Args:
-            u: shape (N, 3), Input array with 3D points.
-            v: shape (N, 3), Input array with 3D points.
-
-        Returns:
-            array: shape (N, 3), Input points transformed.
-            rmsd: The root mean squared deviation.
-        """
-        self.requires_grad = False
-        rmsd = self.fit(u, v, indices=indices)
-        return self.transform(v), rmsd
diff --git a/adaptive_sampling/colvars/utils.py b/adaptive_sampling/colvars/utils.py
index 9965488..297b599 100644
--- a/adaptive_sampling/colvars/utils.py
+++ b/adaptive_sampling/colvars/utils.py
@@ -31,100 +31,6 @@ def read_xyz(xyz_name: str) -> torch.tensor:
     return mol
 
 
-def read_traj(xyz_name: str) -> list:
-    """Read cartesian coordinates of trajectory from file (*.xyz)
-
-    Args:
-        xyz_name (str): file-name of xyz-file
-
-    Returns:
-        [[float,float,float],...]: structure as listimport itertools of atom-data, e.g.: [[0.0,0.0,0.0],[0.74,0.0,0.0]]
-    """
-    if os.path.exists(xyz_name):
-        xyzf = open(xyz_name, "r")
-    else:
-        raise ValueError(" >>> fatal error: xyz file not found")
-
-    traj = []
-    mol = []
-    n = 0
-    for i, line in enumerate(xyzf):
-        words = line.strip().split()
-        if i == 0:
-            n_atoms = int(words[0])
-        elif n == n_atoms:
-            n = 0
-            mol = itertools.chain(*mol)
-            mol = torch.FloatTensor(list(mol)) / BOHR_to_ANGSTROM
-            traj.append(mol)
-            mol = []
-        elif len(words) >= 4:
-            n += 1
-            mol.append([float(words[1]), float(words[2]), float(words[3])])
-
-    return traj
-
-
-def interpolate_coordinates(images: list, n_interpol: int = 20) -> list:
-    """linear interpolation between cartesian coordinates
-
-    Args:
-        images: list of cartesian coordinates of images
-        n_interpol: number of images that are created between two nodes of initial path
-
-    Returns:
-        path: cartesian coordinates of original and interpolated images
-    """
-    path = []
-    for i in range(len(images) - 1):
-        a = images[i]
-        b = images[i + 1]
-        r = b - a
-        d = torch.linalg.norm(r)
-        unit_vec = r / d
-        step = d / n_interpol
-
-        path.append(a)
-        for j in range(1, n_interpol):
-            # linear interpolation
-            p = a + unit_vec * j * step
-            path.append(p)
-
-    path.append(b)
-    return path
-
-
-def find_closest_points(coords: torch.tensor, indices: list = None) -> list:
-
-    coords_tree = coords.detach()
-    coords_tree = torch.reshape(coords_tree, (int(len(coords_tree) / 3), 3)).float()
-
-    if indices != None:
-        coords_tree = coords_tree[indices]
-
-    kdtree = KDTree(coords_tree)
-
-    close_indices = []
-    for xyz in coords_tree:
-        # find three closest atoms (first point is same atom)
-        _, indices = kdtree.query(xyz, 4)
-        close_indices.append(indices[1:])
-
-    return close_indices
-
-
-def get_internal_coords(
-    coords: torch.tensor, index_pairs: list, #indices: list = None
-) -> torch.tensor:
-    coords = torch.reshape(coords, (int(len(coords) / 3), 3)).float()
-    internals = []
-    for pair in index_pairs:
-        dist = torch.linalg.norm(coords[pair[1]] - coords[pair[0]])
-        internals.append(dist)
-    internals = torch.stack(internals)
-    return internals
-
-
 def rmsd(V, W):
     """root-mean-square deviation"""
     diff = V - W
diff --git a/adaptive_sampling/interface/__init__.py b/adaptive_sampling/interface/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/adaptive_sampling/interface/sampling_data.py b/adaptive_sampling/interface/sampling_data.py
index ab7ee7d..0c51ad2 100644
--- a/adaptive_sampling/interface/sampling_data.py
+++ b/adaptive_sampling/interface/sampling_data.py
@@ -22,21 +22,6 @@ class SamplingData:
     dt: float  # MD step size in fs
 
 
-@dataclass
-class SamplingData_QMMM:
-    """The necessary sampling data to perform the enhanced sampling."""
-    mass: np.ndarray  # Masses in atomic units, shape (natoms,)
-    coords: np.ndarray  # Cartesian coordinates in Bohr, shape (3 * natoms,)
-    forces: np.ndarray  # Forces in Hartree/Bohr, shape (3 * natoms,)
-    epot: float  # Potential energy in Hartree
-    temp: float  # Temperature in Kelvin
-    natoms: int  # Number of atoms
-    step: int  # MD step number
-    dt: float  # MD step size in fs
-    mm_epot: float  # energy of mm_system
-    qm_force: np.ndarray  # forces of qm system
-
-
 class MDInterface(Protocol):
     def get_sampling_data(self) -> SamplingData:
         """Define this function for your MD class to provide the
diff --git a/adaptive_sampling/processing_tools/mbar.py b/adaptive_sampling/processing_tools/mbar.py
index 3b1b159..fd904bb 100644
--- a/adaptive_sampling/processing_tools/mbar.py
+++ b/adaptive_sampling/processing_tools/mbar.py
@@ -3,7 +3,7 @@
 from typing import List, Tuple
 from .utils import join_frames
 from ..units import *
-
+import torch
 
 def run_mbar(
     exp_U: np.ndarray,
@@ -12,6 +12,7 @@ def run_mbar(
     conv: float = 1.0e-7,
     conv_errvec: float = None,
     outfreq: int = 100,
+    device = 'cpu'
 ) -> np.ndarray:
     """Self-consistent Multistate Bannett Acceptance Ratio (MBAR)
 
@@ -28,9 +29,20 @@ def run_mbar(
     Returns:
         W: array containing statistical weigths of each frame
     """
+
+    frames_per_traj = torch.from_numpy(frames_per_traj.astype(np.float64))
+    frames_per_traj = frames_per_traj.to(device=device)
+    exp_U = torch.from_numpy(exp_U.astype(np.float64))
+    exp_U = exp_U.to(device=device)
+
     num_trajs = len(frames_per_traj)
-    beta_Ai = np.zeros(shape=(num_trajs,), dtype=float)
-    Ai_overtime = [beta_Ai]
+    beta_Ai = torch.zeros(size=(num_trajs,))
+    beta_Ai = beta_Ai.to(device=device)
+    #Ai_overtime = [beta_Ai]
+    
+    # First denominator with all zero Ai guess
+    denominator = 1.0 / torch.matmul(frames_per_traj * torch.exp(beta_Ai), exp_U)
+    expU_dot_denom = torch.matmul(exp_U, denominator)
 
     print("All ready!\n")
     print("Start of the self-consistent iteration.")
@@ -40,22 +52,24 @@ def run_mbar(
     while True:
         count += 1
 
-        denominator = np.multiply(frames_per_traj * np.exp(beta_Ai), exp_U.T)
-        denominator = 1.0 / np.sum(denominator, axis=1)
-
-        beta_Ai_new = -np.log(np.multiply(exp_U, denominator).sum(axis=1))
-        beta_Ai_new -= beta_Ai_new[0]
-        Ai_overtime.append(beta_Ai_new)
+        beta_Ai_new = -torch.log(expU_dot_denom)
+        beta_Ai_new -= torch.clone(beta_Ai_new[0])
+        #Ai_overtime.append(beta_Ai_new)
 
-        delta_Ai = np.abs(beta_Ai - beta_Ai_new)
-        beta_Ai = np.copy(beta_Ai_new)
+        delta_Ai = torch.abs(beta_Ai - beta_Ai_new)
+        beta_Ai = torch.clone(beta_Ai_new)
+        
+        prefac = frames_per_traj * torch.exp(beta_Ai)
+        denominator = 1.0 / torch.matmul(prefac, exp_U)
+        expU_dot_denom = torch.matmul(exp_U, denominator)
 
-        max_err_vec = np.abs(_error_vec(frames_per_traj, beta_Ai, exp_U)).max()
+        error_v = (frames_per_traj - prefac * expU_dot_denom)
+        max_err_vec = torch.abs(error_v).max()
 
         if count % outfreq == 0 or count == 1:
             print(
                 "Iter %4d:\tConv=%14.10f\tConv_errvec=%14.6f"
-                % (count, np.max(delta_Ai[1:]), max_err_vec)
+                % (count, delta_Ai[1:].max(), max_err_vec)
             )
             sys.stdout.flush()
 
@@ -69,7 +83,7 @@ def run_mbar(
                 "========================================================================"
             )
             print(f"Convergence not reached in {count} iterations!")
-            print("Max error vector:", max_err_vec)
+            print(f"Max error vector: {max_err_vec:14.6f}")
             print(
                 "========================================================================"
             )
@@ -81,7 +95,7 @@ def run_mbar(
                 "========================================================================"
             )
             print(f"Converged after {count} iterations!")
-            print("Max error vector:", max_err_vec)
+            print(f"Max error vector: {max_err_vec:14.6f}")
             print(
                 "========================================================================"
             )
@@ -89,15 +103,17 @@ def run_mbar(
             break
 
     # final values
-    weights = np.multiply(frames_per_traj * np.exp(beta_Ai), exp_U.T)
-    weights = 1.0 / np.sum(weights, axis=1)
+    weights = 1.0 / torch.matmul(frames_per_traj * torch.exp(beta_Ai), exp_U)
+    weights = weights.cpu().numpy()
+
+    torch.cuda.empty_cache()
 
     return weights
 
 def build_boltzmann(
     traj_list: list, 
     meta_f: np.ndarray, 
-    dU_list: list=None, 
+    dU_list: list=None,
     equil_temp: float=300.0
 ) -> Tuple[np.ndarray, np.ndarray]:
     """Build Boltzmann factores for MBAR
@@ -139,7 +155,7 @@ def build_boltzmann(
                 )
             )
 
-    exp_U = np.asarray(exp_U, dtype=float)  # this is a num_trajs x num_frames array
+    exp_U = np.asarray(exp_U, dtype=np.float64)   # this is a num_trajs x num_frames list
     return exp_U, frames_per_traj
 
 def get_windows(
@@ -236,16 +252,3 @@ def deltaf_from_weights(
     p_b = weights[np.where(cv > TS)].sum()
 
     return -RT * np.log(p_b / p_a)
-
-
-def _error_vec(
-    n_frames: np.ndarray, beta_Ai: np.ndarray, exp_U: np.ndarray
-) -> np.ndarray:
-    """error vector for MBAR"""
-    denominator = np.multiply(n_frames * np.exp(beta_Ai), exp_U.T)
-    denominator = 1.0 / np.sum(denominator, axis=1)
-    # sum over all frames
-    error_v = n_frames - n_frames * np.exp(beta_Ai) * (
-        np.multiply(exp_U, denominator).sum(axis=1)
-    )
-    return error_v
diff --git a/adaptive_sampling/processing_tools/thermodynamic_integration.py b/adaptive_sampling/processing_tools/thermodynamic_integration.py
index eb8eb2f..41cb7c1 100644
--- a/adaptive_sampling/processing_tools/thermodynamic_integration.py
+++ b/adaptive_sampling/processing_tools/thermodynamic_integration.py
@@ -1,12 +1,13 @@
 import numpy as np
 from ..units import *
-
+from .utils import *
 
 def integrate(
     mean_force: np.ndarray,
     dx: float,
     equil_temp: float = 300.0,
     method: str = "trapezoid",
+    normalize: bool = True,
 ) -> tuple:
     """numeric integration of thermodynamic force by simpson, trapezoid or rectangle rule
 
@@ -20,7 +21,6 @@ def integrate(
         pmf (np.ndarray): potential of mean force
         rho (np.ndarray): probability density
     """
-    RT = R_in_SI * equil_temp / 1000.0
     data = np.copy(mean_force)
 
     if method == "simpson":
@@ -36,11 +36,15 @@ def integrate(
     for i in range(len(pmf)):
         pmf[i] *= data[0 : i + 1].sum()
 
-    rho = np.exp(-pmf / RT)
-    rho /= rho.sum() * dx
+    if normalize:
+        RT = R_in_SI * equil_temp / 1000.0
+        rho = np.exp(-pmf / RT)
+        rho /= rho.sum() * dx
 
-    pmf = -RT * np.log(rho, out=np.full_like(rho, np.nan), where=(rho != 0))
-    return pmf, rho
+        pmf = -RT * np.log(rho, out=np.full_like(rho, np.nan), where=(rho != 0))
+        return pmf, rho
+    else:
+        return pmf
 
 
 def czar(
@@ -85,3 +89,73 @@ def czar(
     log_hist = np.log(hist, out=np.zeros_like(hist), where=(hist != 0))
 
     return -RT * np.gradient(log_hist, grid_local) + f_corr
+
+def force_error(
+    grid: np.ndarray,
+    xi: np.ndarray,
+    la: np.ndarray,
+    sigma: float,
+    n_bins_tau: int=5,
+    equil_temp: float = 300.0,
+    full_output: bool = False,
+):
+    """Error estimation form the standard deviation in an harmonic bias force
+
+    Args:
+        grid: grid for reaction coordinate
+        xi: trajectory of reaction coordinate
+        la: trajectory of fictitious particle
+        sigma: thermal width of coupling of la to xi
+        n_bins_tau: number of bins for estimation of autocorrelation function
+        equil_temp: Temperature of the simulation
+        full_output: if True return dict with full information
+
+    returns:
+        err_pmf: error estimate for PMF 
+    """
+    RT = R_in_SI * equil_temp / 1000.0
+    k = RT / (sigma * sigma)
+
+    F = -harmonic_force(xi, la, k)
+
+    delta = (xi.max()-xi.min()) / n_bins_tau
+    edges = np.arange(xi.min(), xi.max()+delta/2, delta)
+    corr = []
+    tau  = np.zeros(n_bins_tau)
+    for i in range(n_bins_tau):
+        frames = np.where(np.logical_and(xi>=edges[i], xi<edges[i+1]))
+        corr.append(autocorr(F[frames]))
+        tau[i] = ipce(corr[-1])   
+
+    edges_pmf = np.zeros(len(grid)+1)
+    edges_pmf[:-1] = grid - (grid[1] - grid[0]) / 2.
+    edges_pmf[-1] = edges_pmf[-2] + (grid[1] - grid[0]) / 2.
+    
+    force     = np.zeros_like(grid)   
+    err_force = np.zeros_like(grid)
+    tau_pmf   = np.repeat(tau, len(force)//n_bins_tau)
+    if (len(force)//n_bins_tau) % 2 != 0:
+        tau_pmf = np.append(tau_pmf, tau_pmf[:-1])
+
+    for i in range(len(grid)):
+        frames = np.where(np.logical_and(xi>=edges_pmf[i], xi<edges_pmf[i+1]))
+        force[i] = F[frames].mean()
+        err_force[i] = np.sqrt(tau_pmf[i]/len(F[frames]) * F[frames].var()) 
+    
+    err_pmf = integrate(err_force, grid[1]-grid[0], normalize=False)
+
+    if full_output:
+        return {
+            "tau": tau,
+            "tau_grid": edges[:-1] + delta/2.,
+            "autocorr": corr,
+            "force": force,
+            "err_force": err_force,
+            "err_pmf": err_pmf,
+        }
+    else:
+        return err_pmf
+
+
+
+
diff --git a/adaptive_sampling/processing_tools/thermodynamic_integration_2d.py b/adaptive_sampling/processing_tools/thermodynamic_integration_2d.py
deleted file mode 100644
index dca3c01..0000000
--- a/adaptive_sampling/processing_tools/thermodynamic_integration_2d.py
+++ /dev/null
@@ -1,307 +0,0 @@
-#!/usr/bin/env python
-import time
-import numpy as np
-import matplotlib.pyplot as plt
-import scipy.optimize as opt
-from adaptive_sampling import units
-
-class TI:
-    '''thermodynamic integration with the finite element method to get 2D PMFs
- 
-    args:
-        equil_temp (double): temperature of free energy calculation in K
-        outputfile (string): filename of output
-    '''
-    def __init__(self, 
-        x: np.ndarray, 
-        y: np.ndarray, 
-        fx: np.ndarray, 
-        fy: np.ndarray, 
-        equil_temp: float=300.0, 
-        outputfile: str='free_energy'
-    ):
-        
-        out = open("FEM.out", "w")
-        out.write("\n#######################################################\n")
-        out.write("\tFEM integration of thermodynamic forces \n")
-        out.write("#######################################################\n\n")
-        out.write("Initialize spline functions.\n\n")
-             
-        self.RT      = units.kB_in_atomic * equil_temp * units.atomic_to_kJmol 
-        self.outname = outputfile 
-        
-        # coordinates of bins
-        minxi_1 = x.min()
-        maxxi_1 = x.max()
-        minxi_2 = y.min()
-        maxxi_2 = y.max()
-        
-        dxi_2 = y[1]-y[0]
-        self.x_bins = int((maxxi_2-minxi_2)/dxi_2 + 1)
-        self.y_bins = int(x.shape[0]/self.x_bins)
-        #x = x.reshape(self.y_bins, self.x_bins)
-        #y = y.reshape(self.y_bins, self.x_bins)
-        dxi_1 = x[1,0]-x[0,0]
-
-        # control points
-        self.dy = dxi_1 / 4
-        self.dx = dxi_2 / 4
-        self.minx = minxi_2
-        self.maxx = maxxi_2
-        self.miny = minxi_1
-        self.maxy = maxxi_1
-
-        self.x = np.arange(self.minx, self.maxx+self.dx, self.dx)
-        if len(self.x) == 4*(self.x_bins-1)+2:
-            self.x = self.x[:-1]
-
-        self.y = np.arange(self.miny, self.maxy+self.dy, self.dy)
-        if len(self.y) == 4*(self.y_bins-1)+2:
-            self.y = self.y[:-1]
-
-        self.lenx = len(self.x)
-        self.leny = len(self.y)      
-        self.xx, self.yy = np.meshgrid(self.x,self.y)
-        
-        # coefficient matrix
-        self.alpha = np.full((int(self.x_bins*self.y_bins),), 0.0, dtype=np.float)
-        out.write("%30s:\t%8d\n" % ("Number of coefficients", self.alpha.size))
-        
-        # linear interpolation of gradient to control points
-        origD = [fx.reshape(self.y_bins, self.x_bins), fy.reshape(self.y_bins, self.x_bins)]
-        D = [np.zeros(shape=(self.leny,self.lenx), dtype=np.float), np.zeros(shape=(self.leny,self.lenx), dtype=np.float)]
-        for xy in range(2):
-            D[xy][::4, ::4] = origD[xy]
-            D[xy][::4,2::4] = 0.5*D[xy][::4,:-4:4] + 0.5*D[xy][::4,4::4] 
-            D[xy][::4,1::4] = 0.5*D[xy][::4,2::4]  + 0.5*D[xy][::4,:-4:4]
-            D[xy][::4,3::4] = 0.5*D[xy][::4,2::4]  + 0.5*D[xy][::4,4::4] 
-            D[xy][2::4]     = 0.5*D[xy][:-4:4]     + 0.5*D[xy][4::4]
-            D[xy][1::4]     = 0.5*D[xy][2::4]      + 0.5*D[xy][:-4:4]
-            D[xy][3::4]     = 0.5*D[xy][2::4]      + 0.5*D[xy][4::4] 
-        self.D = np.asarray(D) * units.atomic_to_kJmol
-        self.D = self.D[:,:-1,:-1]
-        out.write("%30s:\t%8d\n" % ("Elements in gradient matrix", self.D.size))
-        
-        # initialize B-spline functions 
-        self.B = []
-        self.gradB = []
-        for center_y in x[:,0]:
-            for center_x in y[0]:
-                self.B.append(self.pyramid(self.x, self.y, dxi_2, dxi_1, center_x, center_y)) 
-                self.gradB.append(self.gradPyr(self.x, self.y, dxi_2, dxi_1, center_x, center_y))
-         
-        self.B = np.asarray(self.B)
-        self.gradB = np.asarray(self.gradB)
-        self.gradB = self.gradB[:,:,:-1,:-1]
-
-        out.write("%30s:\t%8d\n" % ("Elements in gradB matrix", self.gradB.size))
-        out.close()
-
-    #--------------------------------------------------------------------------------------
-    def BFGS(self, maxiter=15000, ftol=1e-10, gtol=1e-5, error_function='rmsd'):
-        '''BFGS minimization of error function
-
-        args:
-            maxiter		(int, 15000, maximum iterations of minimization)
-            ftol        	(double, 1e-10, tolerance for error function)
-            gtol		(double, 1e-5, tolerance for gradient)
-            error_function	(string, rmsd, available is rmsd or power)
-
-        returns:
-            -
-        '''        
-        out = open("FEM.out", "a")
-        if error_function == 'rmsd':
-            out.write("\nerror = RMSD\n")
-            self.error = self.error_rmsd
-        else:
-            out.write("\nerror = diff^2\n")
-            self.error = self.error_power 
- 
-        options={
-            'disp': True,
-            'gtol': gtol,
-            'maxiter': maxiter,
-        }
-        
-        self.it = 0
-        self.err0 = 0
-        err = self.error(self.alpha) 
-        
-        out.write("\nStarting BFGS optimization of coefficents.\n")
-        out.write("--------------------------------------------------------\n")
-        out.write("%6s\t%14s\t%14s\t%14s\n" % ("Iter", "Error [kJ/mol]", "Change Error", "Wall Time [s]")) 
-        out.write("--------------------------------------------------------\n")
-        out.write("%6d\t%14.6f\t%14.6f\t%14.6f\n" % (self.it, err, 0.0, 0.0)) 
-        out.close()
-
-        self.start = time.perf_counter()
-        result = opt.minimize(self.error, self.alpha, method='BFGS', tol=ftol, callback=self.BFGS_progress, options=options)
-        
-        self.alpha = result.x
-        #self.get_F()
-
-    #----------------------------------------------------------------------------------------
-    def error_power(self, alpha):
-        '''error function  
-        '''
-        a_gradB = np.zeros(shape=self.gradB.shape[1:])
-        for ii, a in enumerate(alpha):
-            a_gradB += a*self.gradB[ii]
-        a_gradB_D = a_gradB - self.D
-        return np.power(a_gradB_D, 2).sum()
-
-    #--------------------------------------------------------------------------------------
-    def error_rmsd(self, alpha):
-        '''error function
-        '''
-        a_gradB = np.zeros(shape=self.gradB.shape[1:])
-        for ii, a in enumerate(alpha):
-            a_gradB += a*self.gradB[ii]
-        a_gradB_D = a_gradB - self.D
-        err = np.power(a_gradB_D, 2).sum(axis=0)
-        err = err.mean()
-        return np.sqrt(err)
-
-    #--------------------------------------------------------------------------------------
-    def BFGS_progress(self, alpha):
-         '''callback function to display BFGS progress
-         ''' 
-         out = open("FEM.out", "a")
-         self.it += 1
-         err = self.error(alpha)
-         out.write("%6d\t%14.6f\t%14.6f\t%14.6f\n" % (self.it, err, err-self.err0, time.perf_counter()-self.start))
-         self.err0 = err
-
-         self.alpha = alpha
-         self.get_F()
-         self.start = time.perf_counter()        
-         
-         out.close()
-
-    #----------------------------------------------------------------------------------------
-    def pyramid(self, x, y, dx, dy, cx, cy):
-        '''pyramid function
-        '''
-        return_func = np.zeros(shape=(self.leny,self.lenx), dtype=np.float)
-        for ii, val_y in enumerate(y):
-            for jj, val_x in enumerate(x):
-                if val_y >= (cy - dy) and val_y < cy:
-                    if val_x >= (cx - dx) and val_x < cx:
-                        return_func[ii, jj] = ((val_y - cy)/dy + 1)*((val_x - cx)/dx + 1)
-                    if val_x < (cx + dx) and val_x >= cx:
-                        return_func[ii, jj] = ((val_y - cy)/dy + 1)*((cx - val_x)/dx + 1)
-                if val_y < (cy + dy) and val_y >= cy:
-                    if val_x >= (cx - dx) and val_x < cx:
-                        return_func[ii, jj] = ((cy - val_y)/dy + 1)*((val_x - cx)/dx + 1)
-                    if val_x < (cx + dx) and val_x >= cx:
-                        return_func[ii, jj] = ((cy - val_y)/dy + 1)*((cx - val_x)/dx + 1)
-        #
-        return return_func
-
-    #----------------------------------------------------------------------------------------
-    def gradPyr(self, x, y, dx, dy, cx, cy):
-        '''gradients of pyramid function
-        '''
-        deriv_x = np.zeros(shape=(len(y),len(x)), dtype=np.float)
-        deriv_y = np.zeros(shape=(len(y),len(x)), dtype=np.float)
-        for ii, val_y in enumerate(y):
-            for jj, val_x in enumerate(x):
-                if val_y >= (cy - dy) and val_y < cy:
-                    if val_x >= (cx - dx) and val_x < cx:
-                        deriv_x[ii, jj] = ((val_y - cy)/dy + 1)*(1.0/dx)
-                        deriv_y[ii, jj] = (1.0/dy)*((val_x - cx)/dx + 1)
-                    if val_x < (cx + dx) and val_x >= cx:
-                        deriv_x[ii, jj] = ((val_y - cy)/dy + 1)*(-1.0/dx)
-                        deriv_y[ii, jj] = (1.0/dy)*((cx - val_x)/dx + 1)
-                if val_y < (cy + dy) and val_y >= cy:
-                    if val_x >= (cx - dx) and val_x < cx:
-                        deriv_x[ii, jj] = ((cy - val_y)/dy + 1)*(1.0/dx)
-                        deriv_y[ii, jj] = (-1.0/dy)*((val_x - cx)/dx + 1)
-                    if val_x < (cx + dx) and val_x >= cx:
-                        deriv_x[ii, jj] = ((cy - val_y)/dy + 1)*(-1.0/dx)
-                        deriv_y[ii, jj] = (-1.0/dy)*((cx - val_x)/dx + 1)
-        #
-        return [deriv_x, deriv_y]
-
-    #--------------------------------------------------------------------------------------
-    def get_F(self):
-        '''get free energy from optimized coefficients
-        ''' 
-        F_surface = np.zeros(shape=(self.leny,self.lenx), dtype=np.float64)
-        fitted_grad_x = np.zeros(shape=(self.leny-1,self.lenx-1), dtype=np.float64)
-        fitted_grad_y = np.zeros(shape=(self.leny-1,self.lenx-1), dtype=np.float64)
-
-        for ii, a in enumerate(self.alpha):
-            F_surface     += a*self.B[ii]
-            fitted_grad_x += a*self.gradB[ii][0]
-            fitted_grad_y += a*self.gradB[ii][1]
-    
-        prob_surface = np.exp(-F_surface/self.RT) 
-        prob_surface /= prob_surface.sum()*self.dx*self.dy
-        F_surface = -self.RT*np.log(prob_surface)       
- 
-        estim_err_x = np.abs(fitted_grad_x - self.D[0]) 
-        estim_err_y = np.abs(fitted_grad_y - self.D[1]) 
-         
-        self.write_output(F_surface, prob_surface, estim_err_x, estim_err_y) 
-        self.plot_F(F_surface, estim_err_x, estim_err_y)
-
-    # -----------------------------------------------------------------------------------------------------
-    def write_output(self, F_surface, prob_surface, estim_err_x, estim_err_y):
-        '''write output of free energy calculations
-        '''
-        out = open(f"%s.dat" % (self.outname), "w")
-
-        head = ("Xi1", "Xi1", "error x", "error y", "probability", "free energy [kJ/mol]")
-        out.write("%14s\t%14s\t%14s\t%14s\t%14s\t%14s\n" % head)
-        for i in range(self.leny-1):
-            for j in range(self.lenx-1):
-                row = (self.y[i], self.x[j], estim_err_x[i,j], estim_err_y[i,j], prob_surface[i,j], F_surface[i,j])
-                out.write("%14.6f\t%14.6f\t%14.6f\t%14.6f\t%14.10f\t%14.10f\n" % row)
-
-        out.close()
-	
-    #--------------------------------------------------------------------------------------
-    def plot_F(self, A_surface, estim_err_x, estim_err_y):
-        '''plot free energy surface 
-        '''
-        # Plotting
-        plt.rcParams["figure.figsize"] = [8,14]
-        fig, axs = plt.subplots(3)
-        #
-        im0 = axs[0].imshow(estim_err_x, origin='lower', extent=[self.xx.min(), self.xx.max(), self.yy.min(), self.yy.max()], zorder=1, interpolation='None', cmap='viridis', aspect='auto')
-        cb0 = plt.colorbar(im0, ax=axs[0], fraction=0.044)
-        cb0.outline.set_linewidth(2.5)
-        cb0.ax.set_ylabel('Diff Gradient', fontsize=18)
-        cb0.ax.tick_params(length=4,width=2.5,labelsize=18, pad=10, direction='in')
-        #
-        im1 = axs[1].imshow(estim_err_y, origin='lower', extent=[self.xx.min(), self.xx.max(), self.yy.min(), self.yy.max()], zorder=1, interpolation='None', cmap='viridis', aspect='auto')
-        cb1 = plt.colorbar(im1, ax=axs[1], fraction=0.044)
-        cb1.outline.set_linewidth(2.5)
-        cb1.ax.set_ylabel('Diff Gradient', fontsize=18)
-        cb1.ax.tick_params(length=4,width=2.5,labelsize=18, pad=10, direction='in')
-        #
-        CS = axs[2].contour(self.xx, self.yy, A_surface, colors='black', zorder=3)
-        plt.clabel(CS, CS.levels, inline='true', fontsize=10, fmt="%5.3f")
-        im2 = axs[2].imshow(A_surface, origin='lower', extent=[self.xx.min(), self.xx.max(), self.yy.min(), self.yy.max()], zorder=1, interpolation='None', cmap='viridis', aspect='auto')
-        cb2 = plt.colorbar(im2, ax=axs[2], fraction=0.044)
-        cb2.outline.set_linewidth(2.5)
-        cb2.ax.set_ylabel('Free Energy', fontsize=18)
-        cb2.ax.tick_params(length=4,width=2.5,labelsize=18, pad=10, direction='in')
-        #
-        for ax in axs:
-            ax.set_ylim([self.miny, self.maxy])
-            ax.set_xlim([self.minx, self.maxx])
-            ax.set_xlabel(r'$\xi_1$', fontsize=20)
-            ax.set_ylabel(r'$\xi_2$', fontsize=20)
-            ax.spines['bottom'].set_linewidth('3')
-            ax.spines['top'].set_linewidth('3')
-            ax.spines['left'].set_linewidth('3')
-            ax.spines['right'].set_linewidth('3')
-            ax.tick_params(axis='y',length=6,width=3,labelsize=20, pad=10, direction='in')
-            ax.tick_params(axis='x',length=6,width=3,labelsize=20, pad=10, direction='in')
-        #
-        plt.tight_layout()
-        plt.savefig(f"%s.png" % (self.outname), dpi=400)
-        plt.close()
diff --git a/adaptive_sampling/processing_tools/utils.py b/adaptive_sampling/processing_tools/utils.py
index eea5b56..9f3c42f 100644
--- a/adaptive_sampling/processing_tools/utils.py
+++ b/adaptive_sampling/processing_tools/utils.py
@@ -27,6 +27,58 @@ def boltzmann(u_pot: Union[float, np.ndarray], beta: float) -> Union[float, np.n
     return np.exp(-beta * u_pot, dtype=float)
 
 
+def _next_pow_two(n: int):
+    """Returns the next power of two greater than or equal to `n`"""
+    i = 1
+    while i < n:
+        i = i << 1
+    return i
+
+
+def autocorr(x: np.ndarray) -> np.ndarray:
+    """Estimate the normalized autocorrelation function of a 1-D series
+
+    args:
+        x: The time series of which to calculate the autocorrelation function.
+
+    returns:
+        acf: The autocorrelation as a function of lag time.
+    """
+    x = np.atleast_1d(x)
+    if len(x.shape) != 1:
+        raise ValueError("invalid dimensions for 1D autocorrelation function")
+    n = _next_pow_two(len(x))
+
+    # Compute the FFT and then (from that) the auto-correlation function
+    f = np.fft.fft(x - np.mean(x), n=2 * n)
+    acf = np.fft.ifft(f * np.conjugate(f))[:len(x)].real
+    acf /= acf[0]
+    return acf
+
+
+def ipce(corr_x: np.ndarray):
+    """ The initial positive correlation time estimator for the autocorrelation time, as proposed by Geyer et al.
+
+    args
+        corr_x: autocorrelation function of time series of which to calculate the autocorrelation time.
+
+    returns:
+        tau: Estimate of the autocorrelation time.
+    """
+    lag_max = int(len(corr_x) / 2)
+    i, t = 0, 0.0
+    while i < 0.5 * lag_max:
+        gamma = corr_x[2 * i] + corr_x[2 * i + 1]
+        if gamma < 0.0:
+            #  print("stop at ", 2*i)
+            break
+        else:
+            t += gamma
+        i += 1
+    tau = 2 * t - 1
+    return tau
+
+
 def join_frames(traj_list: List[np.array]) -> Tuple[np.ndarray, float, np.ndarray]:
     """get one array with all frames from list of trajectories
 
diff --git a/adaptive_sampling/sampling_tools/__init__.py b/adaptive_sampling/sampling_tools/__init__.py
index 04a29f6..f32927f 100644
--- a/adaptive_sampling/sampling_tools/__init__.py
+++ b/adaptive_sampling/sampling_tools/__init__.py
@@ -9,7 +9,7 @@
 from .utils import *
 
 __all__ = [
-    "eabf",
+    "abf",
     "eabf",
     "enhanced_sampling",
     "gamd",
diff --git a/adaptive_sampling/sampling_tools/eabf.py b/adaptive_sampling/sampling_tools/eabf.py
index bfbba28..8395f32 100644
--- a/adaptive_sampling/sampling_tools/eabf.py
+++ b/adaptive_sampling/sampling_tools/eabf.py
@@ -323,7 +323,7 @@ def shared_bias(
                         m2,
                         czar_corr,
                     )
-                    self.restart(filename=mw_file)
+                    self.restart(filename=mw_file, restart_ext_sys=False)
                     os.chmod(mw_file + ".npz", 0o444)  # other walkers can access again
                     
                     self.get_pmf()  # get new global pmf   
@@ -470,9 +470,11 @@ def write_restart(self, filename: str = "restart_abf"):
             m2=self.m2_force,
             ext_hist=self.ext_hist,
             czar_corr=self.correction_czar,
+            ext_momenta=self.ext_momenta,
+            ext_coors=self.ext_coords,
         )
 
-    def restart(self, filename: str = "restart_abf"):
+    def restart(self, filename: str = "restart_abf", restart_ext_sys=True):
         """restart from restart file
 
         Args:
@@ -488,6 +490,9 @@ def restart(self, filename: str = "restart_abf"):
         self.m2_force = data["m2"]
         self.ext_hist = data["ext_hist"]
         self.correction_czar = data["czar_corr"]
+        if restart_ext_sys:
+            self.ext_momenta = data["ext_momenta"]
+            self.ext_coords = data["ext_coords"]
 
         if self.verbose:
             print(f" >>> Info: Adaptive sampling restartet from {filename}!")
diff --git a/adaptive_sampling/sampling_tools/gamd.py b/adaptive_sampling/sampling_tools/gamd.py
index 8a52115..b8e3b5e 100644
--- a/adaptive_sampling/sampling_tools/gamd.py
+++ b/adaptive_sampling/sampling_tools/gamd.py
@@ -89,15 +89,7 @@ def step_bias(self, write_output: bool = True, write_traj: bool = True, **kwargs
                 self._calc_E_k0()
 
             # apply boost potential
-            if self.gamd_bound.lower() == "amd":
-                self.gamd_pot = np.square(self.E - epot) / (self.sigma0 + (self.E - epot))
-                bias_force -= (
-                    ((epot - self.E) * (epot - 2.0 * self.sigma0 - self.E)) / np.square(epot - self.sigma0 - self.E)
-                ) * self.gamd_forces
-            else:
-                prefac = self.k0 / (self.pot_max - self.pot_min)
-                self.gamd_pot = 0.5 * prefac * np.power(self.E - epot, 2)
-                bias_force -= prefac * (self.E - epot) * self.gamd_forces
+            bias_force -= self._apply_boost(epot)
 
             if md_state.step < self.gamd_equil_steps:
                 self._update_pot_distribution(epot)
@@ -199,7 +191,7 @@ def _calc_E_k0(self):
 
         elif self.gamd_bound.lower() == "upper":
             ko = (1.0 - self.sigma0 / self.pot_std) * (
-                (self.pot_max - self.pot_min) / (self.pot_max - self.pot_avg)
+                (self.pot_max - self.pot_min) / (self.pot_avg - self.pot_min)
             )
             if 0.0 < ko <= 1.0:
                 self.k0 = ko
@@ -213,6 +205,19 @@ def _calc_E_k0(self):
         else:
             raise ValueError(f" >>> Error: unknown GaMD bound {self.gamd_bound}!")
 
+    def _apply_boost(self, epot):
+        """Apply boost potential to forces"""
+        if self.gamd_bound.lower() == "amd":
+            self.gamd_pot = np.square(self.E - epot) / (self.sigma0 + (self.E - epot))
+            boost_force = (
+                ((epot - self.E) * (epot - 2.0 * self.sigma0 - self.E)) / np.square(epot - self.sigma0 - self.E)
+            ) * self.gamd_forces
+        else:
+            prefac = self.k0 / (self.pot_max - self.pot_min)
+            self.gamd_pot = 0.5 * prefac * np.power(self.E - epot, 2)
+            boost_force = prefac * (self.E - epot) * self.gamd_forces
+        return boost_force
+    
     def write_restart(self, filename: str = "restart_gamd"):
         """write restart file
         TODO: fix me
diff --git a/adaptive_sampling/sampling_tools/gawtmeabf.py b/adaptive_sampling/sampling_tools/gawtmeabf.py
index 88d49ae..41c029d 100644
--- a/adaptive_sampling/sampling_tools/gawtmeabf.py
+++ b/adaptive_sampling/sampling_tools/gawtmeabf.py
@@ -52,7 +52,14 @@ def __init__(self, *args, qm_boost: bool = False, do_wtm: bool = True, **kwargs)
         self.do_wtm = do_wtm
         self.qm_boost = qm_boost
 
-    def step_bias(self, write_output: bool = True, write_traj: bool = True, **kwargs):
+    def step_bias(
+        self, 
+        write_output: bool = True, 
+        write_traj: bool = True, 
+        output_file: str = 'gawtmeabf.out',
+        traj_file: str = 'CV_traj.dat', 
+        restart_file: str = 'restart_gawtmeabf',
+        **kwargs):
 
         md_state = self.the_md.get_sampling_data()
         epot = md_state.epot
@@ -76,9 +83,7 @@ def step_bias(self, write_output: bool = True, write_traj: bool = True, **kwargs
                 self._calc_E_k0()
 
             # apply gamd boost potential
-            prefac = self.k0 / (self.pot_max - self.pot_min)
-            self.gamd_pot = 0.5 * prefac * np.power(self.E - epot, 2)
-            bias_force -= prefac * (self.E - epot) * self.gamd_forces
+            bias_force -= self._apply_boost(epot)
 
             if md_state.step < self.gamd_equil_steps:
                 self._update_pot_distribution(epot)
@@ -161,7 +166,7 @@ def step_bias(self, write_output: bool = True, write_traj: bool = True, **kwargs
         if md_state.step % self.out_freq == 0:
 
             if write_traj:
-                self.write_traj()
+                self.write_traj(filename=traj_file)
 
             if write_output:
                 self.get_pmf()
@@ -172,8 +177,8 @@ def step_bias(self, write_output: bool = True, write_traj: bool = True, **kwargs
                 output[f"metapot"] = self.metapot
                 output[f"GaMD corr"] = self.gamd_corr
 
-                self.write_output(output, filename="gaeabf.out")
-                self.write_restart()
+                self.write_output(output, filename=output_file)
+                self.write_restart(filename=restart_file)
 
         return bias_force
 
@@ -235,6 +240,8 @@ def write_restart(self, filename: str = "restart_gaabf"):
             var=self.var_force,
             m2=self.m2_force,
             ext_hist=self.ext_hist,
+            ext_momenta=self.ext_momenta,
+            ext_coors=self.ext_coords,
             czar_corr=self.correction_czar,
             abf_force=self.abf_forces,
             center=self.center,
@@ -252,7 +259,7 @@ def write_restart(self, filename: str = "restart_gaabf"):
             k0=self.k0,
         )
 
-    def restart(self, filename: str = "restart_gaabf"):
+    def restart(self, filename: str="restart_gaabf", restart_ext_sys: bool=True):
         """restart from restart file
 
         Args:
@@ -283,11 +290,14 @@ def restart(self, filename: str = "restart_gaabf"):
         self.pot_min = data["pot_min"]
         self.pot_max = data["pot_max"]
         self.k0 = data["k0"]
-
+        if restart_ext_sys:
+            self.ext_momenta = data["ext_momenta"]
+            self.ext_coords = data["ext_coords"]
+        
         if self.verbose:
             print(f" >>> Info: Adaptive sampling restartet from {filename}!")
 
-    def write_traj(self):
+    def write_traj(self, filename='CV_traj.dat'):
         """save trajectory for post-processing"""
 
         data = self._write_ext_traj()
@@ -295,7 +305,7 @@ def write_traj(self):
         data["Epot [H]"] = self.epot
         data["T [K]"] = self.temp
 
-        self._write_traj(data)
+        self._write_traj(data, filename=filename)
 
         # reset trajectories to save memory
         self.traj = np.array([self.traj[-1]])
diff --git a/adaptive_sampling/sampling_tools/metadynamics.py b/adaptive_sampling/sampling_tools/metadynamics.py
index 6334c2c..d93aea8 100644
--- a/adaptive_sampling/sampling_tools/metadynamics.py
+++ b/adaptive_sampling/sampling_tools/metadynamics.py
@@ -1,4 +1,4 @@
-import os, time, itertools
+import os, time
 import numpy as np
 from typing import Tuple
 from .enhanced_sampling import EnhancedSampling
@@ -219,6 +219,9 @@ def shared_bias(
                     self.restart(filename=mw_file)
                     os.chmod(mw_file + ".npz", 0o444)  # other walkers can access again
 
+                    # recalculates `self.metapot` and `self.bias` to ensure convergence of WTM potential
+                    self._update_metapot_from_centers() 
+
                     self.metapot_last_sync = np.copy(self.metapot)
                     self.bias_last_sync = np.copy(self.bias)
                     self.len_center_last_sync = len(self.center)
@@ -423,14 +426,14 @@ def _update_wtm(
             new_hist = data["hist"] + hist
             new_bias = data["force"] + bias
             new_metapot = data["metapot"] + metapot
-            new_centers = np.append(data["centers"], center), 
+            new_centers = np.append(data["centers"], center)
         
         self._write_restart(
             filename=filename,
             hist=new_hist,
             force=new_bias,
             metapot=new_metapot,
-            centers=np.asarray(list(itertools.chain(*new_centers))),
+            centers=new_centers,
         )
                     
     def write_restart(self, filename: str = "restart_wtm"):
@@ -466,6 +469,38 @@ def restart(self, filename: str = "restart_wtm"):
         if self.verbose:
             print(f" >>> Info: Adaptive sampling restartet from `{filename}.npz`!")
 
+
+    def _update_metapot_from_centers(self):
+        """recalculate metadynamics potential and bias force from stored centers"""
+        if self.ncoords > 1:
+            # TODO: implement for 2D
+            raise NotImplementedError(
+                " >>> Error: metadynamics currently only supported for 1D coordinates"
+            )            
+
+        self.metapot = np.zeros_like(self.metapot)
+        self.bias    = np.zeros_like(self.bias)
+
+        for _, xi in enumerate(self.center):
+            xi = [xi]
+            if self._check_boundaries(xi):
+                bink = self.get_index(xi)
+                if self.well_tempered:
+                    w = self.hill_height * np.exp(
+                        -self.metapot[bink[1], bink[0]]
+                        / (kB_in_atomic * self.well_tempered_temp)
+                    )
+                else:
+                    w = self.hill_height
+
+                dx = diff(self.grid[0], xi[0], self.cv_type[0])
+                epot = w * np.exp(-(dx * dx) / (2.0 * self.hill_var[0]))
+                self.metapot[0] += epot
+                self.bias[0][0] -= epot * dx / self.hill_var[0]
+
+                self._smooth_boundary_grid(xi, w)
+        
+
     def write_traj(self):
         """save trajectory for post-processing"""
         data = {
diff --git a/adaptive_sampling/sampling_tools/metaeabf.py b/adaptive_sampling/sampling_tools/metaeabf.py
index 240cfc1..6752d17 100644
--- a/adaptive_sampling/sampling_tools/metaeabf.py
+++ b/adaptive_sampling/sampling_tools/metaeabf.py
@@ -1,5 +1,4 @@
-from fileinput import filename
-import os, time, itertools
+import os, time
 import numpy as np
 from .enhanced_sampling import EnhancedSampling
 from .utils import welford_var, combine_welford_stats, diff
@@ -285,9 +284,12 @@ def shared_bias(
                         delta_metapot,
                         new_center,
                     )
-                    self.restart(filename=mw_file)
+                    self.restart(filename=mw_file, restart_ext_sys=False)
                     os.chmod(mw_file + ".npz", 0o444)  # other walkers can access again
 
+                    # recalculates `self.metapot` and `self.bias` to ensure convergence of WTM potential
+                    self._update_metapot_from_centers() 
+                    
                     self.get_pmf()  # get new global pmf
                     self.metapot_last_sync = np.copy(self.metapot)
                     self.bias_last_sync = np.copy(self.bias)
@@ -323,7 +325,7 @@ def _update_wtmeabf(
             new_m2 = np.zeros_like(self.m2_force).flatten()
             new_abf_forces = np.zeros_like(self.abf_forces).flatten()
             new_ext_hist = np.zeros_like(self.ext_hist).flatten()
-            new_centers = np.append(data["center"], center), 
+            new_centers = np.append(data["center"], center)
             
             for i in range(len(hist.flatten())):
                 (
@@ -348,7 +350,7 @@ def _update_wtmeabf(
             ext_hist=new_ext_hist.reshape(self.histogram.shape),
             czar_corr=new_czar_corr,
             abf_force=new_abf_forces.reshape(self.abf_forces.shape),
-            center=np.asarray(list(itertools.chain(*new_centers))),
+            center=new_centers,
             metapot=new_metapot,
         )                           
 
@@ -368,9 +370,11 @@ def write_restart(self, filename: str = "restart_wtmeabf"):
             abf_force=self.abf_forces,
             center=self.center,
             metapot=self.metapot,
+            ext_momenta=self.ext_momenta,
+            ext_coors=self.ext_coords,
         )
 
-    def restart(self, filename: str = "restart_wtmeabf"):
+    def restart(self, filename: str = "restart_wtmeabf", restart_ext_sys=True):
         """restart from restart file
 
         Args:
@@ -389,6 +393,9 @@ def restart(self, filename: str = "restart_wtmeabf"):
         self.abf_forces = data["abf_force"]
         self.center = data["center"].tolist()
         self.metapot = data["metapot"]
+        if restart_ext_sys:
+            self.ext_momenta = data["ext_momenta"]
+            self.ext_coords = data["ext_coords"]
 
         if self.verbose:
             print(f" >>> Info: Adaptive sampling restartet from {filename}!")
diff --git a/examples/JCP_2022_mbar_and_eabf/double_well_potential/analysis.ipynb b/examples/JCP_2022_mbar_and_eabf/double_well_potential/analysis.ipynb
index 428a6df..47d1519 100644
--- a/examples/JCP_2022_mbar_and_eabf/double_well_potential/analysis.ipynb
+++ b/examples/JCP_2022_mbar_and_eabf/double_well_potential/analysis.ipynb
@@ -50,7 +50,8 @@
    "source": [
     "# run MBAR and compute free energy profile and probability density from statistical weights\n",
     "traj_list, indices, meta_f = mbar.get_windows(grid, cv, la, ext_sigma, equil_temp=300.0)\n",
-    "weights = mbar.run_mbar(traj_list, meta_f, conv=1.0e-4, conv_errvec=None, outfreq=100, equil_temp=300.0)"
+    "exp_U, frames_per_traj = mbar.build_boltzmann(traj_list, meta_f, equil_temp=300.0)\n",
+    "weights = mbar.run_mbar(exp_U, frames_per_traj, max_iter=1000, conv=1.0e-4, conv_errvec=1.0, outfreq=100, device='cpu')"
    ]
   },
   {
diff --git a/examples/JCP_2022_mbar_and_eabf/skewed_potential/analysis.ipynb b/examples/JCP_2022_mbar_and_eabf/skewed_potential/analysis.ipynb
index b67ce5e..a04ae14 100644
--- a/examples/JCP_2022_mbar_and_eabf/skewed_potential/analysis.ipynb
+++ b/examples/JCP_2022_mbar_and_eabf/skewed_potential/analysis.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "id": "f18eb4cb",
    "metadata": {},
    "outputs": [],
@@ -25,7 +25,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": null,
    "id": "ec7c7d81",
    "metadata": {},
    "outputs": [],
@@ -35,7 +35,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": null,
    "id": "c5c713d6",
    "metadata": {},
    "outputs": [],
@@ -72,7 +72,8 @@
    "source": [
     "# run MBAR and compute free energy profile and probability density from statistical weights\n",
     "traj_list, indices, meta_f = mbar.get_windows(grid, cv, la, ext_sigma, equil_temp=300.0)\n",
-    "weights = mbar.run_mbar(traj_list, meta_f, conv=1.0e-4, conv_errvec=None, outfreq=100, equil_temp=300.0)"
+    "exp_U, frames_per_traj = mbar.build_boltzmann(traj_list, meta_f, equil_temp=300.0)\n",
+    "weights = mbar.run_mbar(exp_U, frames_per_traj, max_iter=10000, conv=1.0e-4, conv_errvec=1.0, outfreq=100, device='cpu')"
    ]
   },
   {
@@ -135,7 +136,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": null,
    "id": "cdf887bf",
    "metadata": {},
    "outputs": [],
@@ -145,7 +146,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": null,
    "id": "46f5b5fb",
    "metadata": {},
    "outputs": [],
@@ -157,33 +158,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": null,
    "id": "c8f82862",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7fc7b8790730>]"
-      ]
-     },
-     "execution_count": 15,
-     "metadata": {},
-     "output_type": "execute_result"
-    },
-    {
-     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEGCAYAAACAd+UpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABVr0lEQVR4nO2dd5wURdrHf8/MbGATu8CC5CULgiiSTJgwYnzNOWc9PT09zNnz1FPPU+8OxTNyHCbABCIYMIFEJYOwRMksLOzu7IR6/5iwE7p7OlR198zW9/NRdrqrq6q7q55+6qmnniLGGCQSiUSSm3icroBEIpFIxCGFvEQikeQwUshLJBJJDiOFvEQikeQwUshLJBJJDiOFvEQikeQwPh6ZEFE5gNcA9AfAAFwFYDmA/wGoAlAN4DzG2C6tfNq0acOqqqp4VEkikUiaDXPnzt3OGKtUOkc8/OSJ6E0AMxljrxFRPoAiAPcC2MkYe4qIRgOoYIz9WSufwYMHszlz5liuj0QikTQniGguY2yw0jnL5hoiKgMwAsBYAGCMNTLGagCcAeDNaLI3AZxptSyJRCKRGIOHTb47gG0A/kNE84noNSIqBtCOMfY7AET/bcuhLIlEIpEYgIeQ9wEYBOCfjLGDAewDMFrvxUR0HRHNIaI527Zt41AdiUQikcTgIeQ3ANjAGJsV/f0+IkJ/CxG1B4Dov1uVLmaMjWGMDWaMDa6sVJw3kEgkEolJLAt5xthmAOuJqE/00HEAlgCYDODy6LHLAUyyWpZEIpFIjMHFhRLArQDejXrWrAZwJSIfkAlEdDWAdQDO5VSWRCKRSHTCRcgzxhYAUHLfOY5H/hKJRCIxh1zxKpFIJCbZVuvHlEWbna6GJs1ayO+uD+DxT5agMRh2uioSiSRKOMwQCGVHn7z89dm44Z252OcPpp0LhZnicbtp1kL+manL8Np3azBx/kbLee3c14gBD03F/HWakRskOUY4zBAOy93VeHL5f2aj132f21rm3LW7MPSJL7GnIWDouvU76wAAIcZQ2xDA77vr4+cenrwYBzw01fEPVrMW8jENPqwQ2qEhEMLWPQ2685q1egdq/UGc9coP3OqXiZdmrMRv2/baVp4T/LKhJt6R3EiP+z7DSX//1ulqZGTHXj8++WWT09VIYt2OOlSN/hRVoz9NqtvMldsBAD/+tiN+jDGGjTX1aXlo8d6c9aga/Sn26tCmn5+2Altr/Vi4vkZ3/peOnYXahLxHvfgdDv3LjKby564HAARDzioBzVrIx2Q7Ufq5S8fOwtAnp5vKV6/5pyEQQsikFrhq6148+8UKHPe3b0xdz5uGQAhPfrYUdY18h6env/Q9jnz6K6558oQxYMWWpg/tWz9W47WZq/HEp0uwYkstACAQCuPlr1ahIRByqpq4+s05uGXcfOzY63esDqn89+d18b8nL0j/AH23qmlx5Ds/rcXhT83Arxt2687/X9/8BgDYvLse/mBI9dpAKBx/N0ZCecU+RrH6rUtRRtyyfXbzFvLRfwnpUv7navNmF6WPhhL7PzAFd05YYKqMxA/J0Ce+xJRFv2e8pmr0p+j7wBRT5WXijR+qMebb1Tjl7zMzpg2EwoZNHFWjP8XaHfvMVo8bjcEwLh07C4s27saa7en1eXDSYjz+6VK8OnMNLh0bWR84/uf1eGbqcrzy1SohdWoIhPDR/A3Yta8RVaM/xaQF6ebHmBZsVqkQDQMwp3pn0ujZH2hqJ7PW7AQAzFimuKYyjYZACLFbZQy44e25OO2l7/CXz5aiavSnuG38fHy5ZAv8wRDO//ePmLM20t+rdbSxQCicZJYBgKenLFdNr1ceiIKXn3xWwRjDIx8vwd6GqNZJEdtqYyiMwjwvh/wj/w578kvcckxPXHpoVVqa76JawMQFm/DCBQdbKm9rrR83vDMPd53YBxcM6YzWJQVpaXbXR2yN9YK0yd+jQqR6Rx1WbKnFNW/OwdtXD0XX1sVpaWP21k9uPQLFBT54idCldVHGMmat2amYn530vj9S95krv8uYNtYO6qJD+rpGMc/+lnHz8OXSrThlwH4AgNvGL0CvtqXo16EsrS7b9vrRsigPBT7r7dwqK6MjHQCYtmQLpi3ZknT+te/W4LXv1uAPx/bEV1Hh/vyXK3DbyF4Z894/QZlhAL5aHhkV/Pvb1QCASQs2YdKCTbj6iG6Yt64mnvbBSYtxmUJ/jfHXKcvwz69/y1h+rFw3kLOafGMwjFvGzVPUtl6asQpv/FCNKYsjrk/TlmzB6A9/wf4PTEnTdH74bXva9QAw/MnpeHjyYlXzxPqdddiyx48HJi1OO1fbEMAlY2cpXGWNZ6Yux93v/6J4buAjX3AvL4Y/GMKbP66N/z7h+W+xbmcdjnrm6/ixUS/OxF8+W5p03an/+A7HPPs1RjyTbo5hjOGGt+cmHXPzBOcRf52BD+dtSDoW0+DGRAXLrxsj5oIBD0/FPR8qvyczfLk0IgA/+7XJle+UF2dia21EK35g4iJsj5ppRr34HYabNEPyJlbvTLw4YxX2WfhAaplNPv818wg4hj8Y0i3gN+9u0DTbvvzVKlzz5s+6y7ZCzgr5OWt34pNffscxz36dJoj/Nm1F0u9pS7ZgwpxIB122eU/Suc9/VfaB3bynAW/8UI1+D07FkU/PwJaEYea6nftU7ciPfrwEAx62LnCZip6wj7NNXA9fLN6SMc3iTXviWpQewgzxj3CMlwSZO/Si9ZHZsKsed0xYmHRsy56IYN2xrzGe5rNff0dtQxD/nb1eNa/3526Ia65WGPrEdAx+fBre/mlt0vFddQGs2pq7E/ZGTFKbdut3rvjTe/o+zOEww/C/NH1IZ0dNTQAwacFGTFn0O56Zulz3R84qOSnkJ87fiItebdKUx81ap5E6mee+WJExTerih/U765ME0uJNe1IvifP692vSjlWN/lR3/WLsrjPm6pUIz3UBa7bvw63/na96flOKR8RqDW8gxhien7YCW/Y0QGkzmw27jHlX8OaeD3+1dP3Gmnrc9O68+O9FG3fjjJe/R32Klvqn9xbiyjf4aHnb9zYqHh/5nDsm7EWwp95839BixtLMygzQNGKLkdjfbhu/ADe8My/1EqHkpJC//X8LTF9LKbMkqVoQANzwzty0Y9tqm7wWAiouU2qmHzNc9JqyuSesQ37znAjKtNjjkrGzMGt1kyvcsSreQKe/9B263fMZ/j59JW4fv4BfBTnyvznq2rcaWh/jxz9dgoXraxxbW5H4XnIZtVGvKJ79Qn0S1glycuK1HXbiD76PcI73W0wKHYbd4eeTzh/rmYfX85/FFlaOD0IjsJG1wVrWDsd55uFndhsA4HTPD3gx/yWc6n88Lf9OtA1HexZgkGcl3gsdhUu80zBm+6moLnwQNzf+AYSBmJF/B7p7NuOghn/Hr7tx7Nc41zsLjcyHzaw1/lfwGPaxAtwduB7AKEP3eCD9ho/yH8Rn4WEY5FmJjrQD44LH4IO1IwAcmpZ+P+zAf/KfxnLWGYydbKgsLTxEGJv3DAZ41qAt1aSdf2bnefj3a13wcf77mBfuhXoU4gbfx/gkNBw7WCku903D84GzseX3ClR5CuBDCBsCkfrd6J2M/p6IiWeUdzY+Dw3Bxg0D0LFTF271N8KhnsX4b/4TWBNuh/9rfAT/y38MjcjDblaMjawNDvGsQA/P73gocDkeyXsTjwcuxj/G1+Efea/hNO9P+DHUDx+Gj0CYefBTuC9KUISzPDMRDA+NFBAKAKumowdtRAPL112vAjRiTN5zuDrwJ5zimY0rfVMwkH5DEB5c1Hgf/uj7AG1oN5ayLjjWswAMQEuqw6TdiwC0FvKs9HCSZzb+lf8CPgkNw+TQYfizbzxKqT7ejjaxVvgx3A8Lwz3waF5kk7m/B8+Cnr5SXXhRNP3/YdacIpzt+RZ35f0P61klhnhW4LXgybjG9zmuaLwLV3s/x4/hfjjZOxtvhk5Uzb91YBMWF/4Rc8O9cIhnJe4JXI2/5I2Nn/8kNBzfhw/A/N2nYkb+HQjDg22sHHsxMZ7mZu9E3JU3Aaf6H0cdCk09N6Nw2eOVF9z2eH24ZdLPOd2ux+DLn1Y9n8g3pafiqDvfTU7zcPLwK/hQBXykrjL/MPQlHDb7lrTrZz04DMM8y1TqrN//N5Je/R4U80pIv/X6RWjbvrOx8lT4bdlC9Bg/gkteMT5rcSqOv+YJ5P1jYNq5NeiIbg8v4VqebrSeuQlqUYRS1OHzQ8fh5BNHYeeke9Fq/ssJ5elrE7se6ogK2ostrBztFD60aswtG4lD7vjAYK05YvZ5Znguu3btQsXfq8zlrZW/zvp+lnciTglMjf/++ag3MeSYM5XzMNrvVRC6x2s2UL6vWnfailDm3am0BDwAlO1VnmA8mFbqrodIAn79k02Z8AX5r0btElynalPqBushKNxCKSLPbuHqyD2tWZnuiaWHCorMcxgR8ABQXmOuPPfjrOLaPZw875bn36mS0h6ahZAPMf1G6IqguC0IHV4T0QRPo7yHfxMiMLDm0TQBAId3bwUAYDa3kLY+94aLsAQ5vwYgEdKYKAvt5TdPp0az6Ek1BmbbKwP6/WbVcY04F49HxLSOe0yIdpDvi3RDu++aCXl3klSIqTsn1G9Pd+zgTbMQ8kbgsd6GkfJjJdcIL34fIfLyFxQE+z0iHCU+L2bvPQcpz9by7MPZtpMa8HBHrYbrr4qs4EmzEPKGhsHMug85UzGHuEfI84M8IgQFuSa4kx3E2gUzYFbkQYMncyiJbMTpWDEeJMuQUEjDzdgj3rQkhbwAKOWxjnp2SvQ4H6oVQjUYg6cm7y77ZzYS83Cz+7u2I6+9zSXag9MKwv7h5NAHyzfvwYSf1yuHBbdBk28WRjm7Q56kflQ+3Xs+AD6uUgCwYeM6VHHLzRoeEZNcxOD0kNtOiDHUNQYR1LOSjSN2Kz/NlfU76/DMB5GQCNWprvHSXMOH3/e4I4a2h9QFl5F4G2U7rAW3UjMnuYrmI+PBwDDosWnYoRKCQFy5WdAOTOC25l2BvWkmnDhSyPPByPCNh91cbeJV8xqnx5gmIQ//HuWyPmoDDC0D29GJxLnvKtF3H/9IqK7AZX3pz3njMTn/frRHehgJJm3yfAgbuE0PDxXSzNfZwIQvWfz6t+AQMz+GKIHsrm4qFmJhzCq8BYM89kbZJMYwT+5JbAv9PdX4ouBuhTNSk89KhAsoi+NRn5fjxKsAKZ+LXkharNrCb77GCPUsD//3yg/Ytc9eM1FzpZTSXSmZDQu3moWQt32RiQn9lnFw3dSN24yWCrhsxC2U43/7iyPlbmdlOICq0RB0bu/ZZo+0yfMhbEDoqk6QGMDUhFZzkmo6aE5Po5Kc0eS7ezbj04J7QX71/Q+yk+xpPSQgLEgq3EogIi8RzSeiT6K/WxHRNCJaGf23gldZJmqnOyWP5lH+e+b9P9PKrdff0a0q4mu21WZOpBMR3yYSlbFEEQq5w/usOWJHjCaeJdwGIHETz9EApjPGegGYHv3tengYMjptM77rDgsZiQxprZbhde73qshWb6PsxP3mu1xl91p+e/2qwUXIE1EnRCLtv5Zw+AwAb0b/fhPAmTzKMoPd/sAhZsKFMq9UQE0iJO5aBfBeHCaFcbZjxuXXzWSTflCzS3wYYl5v9wUAdwNJBu12jLHfASD6b1tOZZlCaxPmRHjY5ENmHqvAydDUzcmzoRNkQRVzhpo6+zd/F4mItvNztRhhzASsM0nFspAnolMBbGWMpW98qu/664hoDhHN2bZNzGKQi33TMXlKJH4MguLtj2aE/OfvPJ85UZTUfWgzpm9MjnXjdgFKzSsGpeMs3ODMxK8oRJj6ausFuZlmiQvl4QBOJ6JqAOMBHEtE7wDYQkTtASD671alixljYxhjgxljgysrKzlUR5l2aycBAPw7jW/GbJQgjL+4sza/YCC1MSFfsn1h0u/FG93tTZG4u71EPLn2tPM4rgOJUb6dw7akSmSDCyVj7B7GWCfGWBWACwDMYIxdAmAygMujyS4HMMlqWVaIzWJvrdXW5Hk0DyMrbG0hRbPptPUrhyqij4Zg2P3DjVwix2zyIgwgvrAoTT4LhLwGTwE4nohWAjg++tsx9E4uaQUR04twv2eD5hqWoqt1CPPY/Uosti4Oa+bkaqAynnhEzZllW6hhxtjXAL6O/r0DwHE887eCkYYcDjS6TRdPwqjJkaVMOrtdSW5uYQ2cxu7NSkQjovUI84vIck3eVRhxEws3Wt2UQyy//fyZofRpWjHPiSlRrjp2y/ntKwG/wqYOzYCsCD2dyH9GAZt/tbVIUZq81WCDemg+Qt7IkNTlNsoRtZ8buyBFyGddp7aDlwYD4853uhYOkWXtYe132DrhdtXT3upvuRfJJTqtAiGf+C0Y3S3NuKL/VutqFB2BspZUZZvMat/hMBBIXpkrwl3NMRfKtcbDUeQC2bBuIpXte9UdKLZ8/y738kTpRcFC8dFemo2QN6K9btm8UWBNeGDYKM+n2I9vBZ5oxycvDQZ5VqF8zCDh5UgiBLJQyGuNzHfU8feEMbo2RS9zq3fpXqhplmYk5PXfqtvjphhtE6n3Y/ru5r9j9kqJi9lvp6l1jA6jLnQDQTGjSxG8PWstXp25WkjeMZqNkDcUiZI1xddeEO4hojKWMOry5vaPlsRZSv2bnK6CYbRG5jv28V/VLsyFEsCKLWIn/JuNkDeiydc2BOJ/F3LcKs8pUi3cpQXZf0+m2LYC2GvvPqrZgB0eHrzR0ltE+P2LEvG3+z7Ezp3pe7/yhKufvJspbdTfufck2PR8Lmz/hoeOKTb5/Lxm89qTeXkIkF8C3Js858JY814OZMfGFdyx2UNMVHHneL/F0D07APwgpgA0I01+8Nb3Aej8ImvI0JDgSRI9GBVJ21NCOexsUWW5Dvd+9Ct+WLU9+0xBjc3TF14LFzRpE9gr5D0Cy+tSt1hY3kAzEvJGCGt4o/z729/MZeqgMNyV4m3AoyrjZq3DRa+5f/MRSWaycRyjVWcRk6TEIdyJUzTTcbs2SStEUyTihh0mNUGOQt5oI25dnG/peiWqCy8CAOysdjQkER+ybTTCnewT8lp1HkBruJdmQ9h3YTQ7Tb7kq/szptHyWz1863iTJTsnSCqK8pJ+81xqVP79o9zycgo9Mj7rzFIGyMY70wq3U+XZwr08UX7yALAg3F1Y3kAzFPLlG2ZkTOP3lqiea+uvNlcwRyHRmrQ34m4MhJJ+v/FDtbC6BJpJ7PcclvFZKeTtHn00BkOZE5lkXri3sLyBZijkM/FbuD0a8soF5GxfV6KU0MLpWii/ugSzc9YumZ1iF6O4HV8WetfYPY+wcL248OGio65m39s1SdhAowhpPPNAyKTm6qAqeIp3dtJvI6aHG9+Ziw/nbchpcwV++zJjkhy+e1S1KXa6CoaxW8hLm3wWYCisQbhJkO/0tEo6t6mm3mwNTF5noqSUogZ7lpvO6/NFm3HHhIXw57RZRkcPzuGPXDZ619htrinO4gWEzUbIN3jLdKcNJ3Ro8iQ7IBXAZPAjG3c6Sp1YTQ2TumLLXuz1B3XlNZiWoR12AoGEj1uOCbzNuxsyppE7VbkLu/c5ad+yUFje0lzDiXntztGdNtFPPj9lyetp3p/MVcBOIV+3M+l3e0r+fY73W2zepT15G+P9gkcxreBusGBTqIdE00026oCprN9VlzFNbn3WUsjK/QXE1nnzom+Bh1tiy5JI+OksdpNvPkI+TPqHW5UlBU0/eGmtdmqC4cyeAHm7q3VnV0Z1GPvqC/HfOabI6xNyOXfT2Y1oE9Oanz4CAKye9XG0PHHvX2ryvCDC8o36AgF5NZ7K4nBXc+Xbaq7h3wFuqX0hIf9cI/PzyumJZ6Osnw007FE9PfOTtzDrq0k2VkggOfDam42QJwB9XtW56ICp/sA2Vm6uAjYKibU7dOxRa6E+ifbpHOgDSezQ2HEoVzG0kXfjPmDs8cD/LlFNcuScWzHsm8s41EydNEVm86/AhMuAkL65pkw05R5t6wL7r2hjWbMR8m3q9S11JrCkiVdehMIWNfnaLcDDLYF5b/GpkAXxnCjks9Gam07TXagFoMvpiVcDNnkWjDgehDbOF1UbfaTW+YNrgCWTgB0rnamPBS7yTheaf7MR8npNGKlCnpe97NVvV1nLYGc0MNqCcRmTGtLMUli1tRb/+kY7CBtLEIS5oMnreVy5ba7Rd2+NjQFs/etBAID6Rj4aMy9iLr68o8Ta4V7qI7EKRLMR8juL9JlqGCipQ/fbl7yQyANzL2TVZmsr5mJ10rMYK+QryJhGrene8spErJg6Bn6NZdxhHRO7WYvag8lhIa9XkNXV1qAd1YitjG6S67x5d8TFd+1OHaZKE4j8yH8VGigsb6AZCXlDJLxQH4IRU0mUQ71LTWWpZ8Uc86u7Na7YEjm3ZrueRmzeW+QN9gCey/8X6hsasXqbWsRNh10oA5n92iUG0Cm/tEJw203qh6mEid0noHT7PGF5H+NdKCxvoDkJeZ3SSMkm729oEqweZk6L1SPkt+7YpXouEI21oGc4albpYIyhEjUAgEfGvo+3XrhHOV2YvznLECEnJkdzV5NPYta/I3M/jelrB9xssmrN1PuOYRjD0t/3xP8GgLZrJvLL32aaj5A3QKrCkrgdoFkbna6PAye1WJ/GxcAYw1+nLMPiTRFT0rNfNIU/eH7XLXg4T3mSl7Eme6wjS+K5a5RN3cCj4hrIciEQmwqJm2Kz716I/FGX7m6cv2xi/O8S1AHbnZvkFNnuQmGGU4PTUo5mr4uBZSFPRJ2J6CsiWkpEi4notujxVkQ0jYhWRv+tsF5dC/XUqYkxUNrCB5YgoM1qrifUmI1DH6tXtHw9lhidefqDYfz4zVRc+kokQFfhj8/Dq2Npn9MCjwUFavLhgOJhkYth1GgI2DX30dSoauojCs2ehvSyS6f/OfnAS4OF1koT1Y5g/T2xUACVFFF8Yh9AI7Gv3AaPmgcB3MkY6wtgOICbiagfgNEApjPGegGYHv2dFWiF5jX7Pe9R/6uOctXPkYGhcr7Waq5YfmBAoB4TCx7Ey56/AQAuxqe68mdw1lwTMhsJVAfqH1H777OmTvmDEw4zzF/HzzyReGex/QHqecRPf+ds63moIHQEmbQzXOSfnR2OEleeYCwLecbY74yxedG/awEsBdARwBkA3owmexPAmVbLsoYBX+AUcwBZ9XEHUFuv3GGTytWRj57GvV9ZZu8axhiIReo0gIzFU2ccnocVQmH73fecHr0kMva7NTjrlR/w/artwsrgYn5flTmEs3nU+oF14a906zWVQy3n6xRcxyBEVAXgYACzALRjjP0ORD4EANqqXHMdEc0hojnbtm3jWZ0U9Lfa1O3/KEHoe0xGKkrdyMMocXONjvvQW0My2SEShbwTNnnGWZPXcwdOiHi1UcXKzbswyvMTNuoIrKaHxHuLFclcHrRMvd1Zf1PJc1qR/Nz+PLTgJuSJqATABwBuZ4ypB7ZIgTE2hjE2mDE2uLKykld1TKPceKwLlb6e9ZnL1mifvB0brJhZ3vmpmks+puHop3/BmB/x8S+bdKR0jyY/cttbeDn/RXTYknkrS30ktvmoUOPc4K5/ew6qRuszBzqN8qitmQt5IspDRMC/yxj7MHp4CxG1j55vD2Arj7LMo/8lpb1jV7mOcYyYGE1n9O5aLP3A4BV8SfTuscpBa9/AIbQs/ltVP3TYRJVIy0CkKxX4a7jkp6TJ8/6mbVsyE1d4p3DLrygs1i8+Rkzpc2LinRe+zEm0ocg25mMBLGWMPZdwajKAywE8Ff03e8LSpdjkwzZ1cC3Th8hGZjTnYXVfi6iGbhhHTX50njWvp9xAvJb6YcHD0b/+wSW/A+tncclHmewV6Erw0OQPB3ApgGOJaEH0v1MQEe7HE9FKAMdHf7ueAgqkKcJuWsbPqzuyMEv4cGTZUJSTn/wPvylNXCo/i+rPnudSphGcfCvZuSUgH5JHbTHzVfY+D8uaPGPsO6i3x+Os5m83nWh7mndNnie7vux67Kl/em8B3v5Tb8tlpdrkA8yLPBL8UeQ0shoz9l84LD/5GKlMsPVd9pJrvoUiqxHL21u7EdjlASqqBJbmTpJWvJs0abqJ7PXwF0iqTZ5MhjLgSaxKvBSKYcE5aXnrFR/1ti3SUYGTC+Ub+c+kHXvzh2pMmJM+Se7IBLNNMIVf+713KvB3sYGz3IqbwzeYoRkJeSMvLsWF0i6bvI1t68689xXK01eBRIGXOqy34xaYipC/bfz8uAfHvR/9ij/+b4HhvC+edx4qJ12cdtwlSrwgsvfuvlq2FZ/+8jvXPJPMNbEVr1ks+C2ba7KF2Eo+PaS+z137GmBHTAatyVVDjUx32uR0Pg6uogweAII1fRWb/KQFTa6QLee8hPbUAJz/pmJaNSppj2JUQLPrIyyRSfaK8N0mZJVt4so3fgYAjCrkl2fi7Td1pez9EDYbTX5Btf6FVqneNL9tMbeEfGmLQQav0NOQjDW2yW89p34ybm+M5FkG47G4HTFj6Jh4/XPeeNzqmyi+LgKhRm03QV5PPnE0lm2ijBBOWmjI41uspFCV7VpkPWOHyBkh/5fPl+Lmd9VjPpeTEb/a5JdcunetqTqFOT5eY8amptSnr35ER3rC9gWf6c7fcUGg4u10uud7POEbqyuLqU+Ii6vCC2pU/ujyfv5KNvls4Zv8P2JxwdUJR6zXP9E8G1NivtmtuGA/K8gZc82/v4nEX3lZ5fyVvqm680qNJ19eay6kKjMo5PUEKNPXwY039DYTL9SdlgxP1fJFbb/VF/PV3n46JwZExlXhw+76AOYs2oyT+u+XdDz+9HnZiXnKdZtt1108ySN0HuEHlMymCzbVZ61KnKXVTuci73T8wfth5oQ6SB2uxQJ5GSVsMDyp/TqU2RKbriug5Gdji/nGResWRHLXewtxwztzsX2v2E1SEgWjWRG5qaYe/mAI7Gd9IylR8DDX1Oxr2nksZspK/dBmEzmjyT+ZF2tc1htZqjKyvbQv+pjKyWCX0dBChCxOMekD7LS5hpSEfA5uCbhv9w4c5dmOYCh5uUn8+XOaeE1s72Y+0o3BMI596nMcP6ALbg/8jB5cauUcK8f/GV2jf6/fGQkCV1maObKrW8kZTZ4rKVI+4GmBXzbUGHajMqrJc8PgilCjH5CiPIebjYKQ979zvgMVEcu/8v6GN/P/Cmqosa1MPdtLphIM+LGs8EocsvSv+H5V+o5SSohySYxtk2mF3g1Nez/8vtse5WHq4s3C8pZCXoFwijazeP12PP/KS/hw7gZD+Ri3yfNxoSxYoy86odl+1hBQX4xkx3J4pQBlBWu/jv+dGiraKuEG/ZP2fpbHrdxuFOn4Db9OBB5uCf/O5EVarjHJhyK7SV3h+wKX+VK3zVNm66Zq1XOrC/bPeP0+pqxZ/+XzJbrK18KjMOeUP0+sGeq+t2dgT4M5s3AmpJBXIFWgHlo7Bf/JfwbFi982lo/B4bS2n7yBjIJitQ+P0x4YGRan/fCbPm1SL4vGP6A77auhU7iWDQD1P0YEzKKfv006nsUhzuHV2Ix9c9mBpvPdWmM9xr5S6x4QyLyzmxXmFN6IUL2YyJpSyCuQKmw7RaMkF9X9jgDzOlGleI/m26+T/eT1oiXkbRE8GSZepy/lO/St21erO219h0O5lg0AQZXHLUKTb0367zV+/ZunG76GLAaZKyY/MDN9DUgpWRfyySu6bSRYLyRbKeSVSHmzS0OdAQALN+4xODFlTOIpTijG66S/XL0dKOaKaLQh9/WsUzy+NtwW6w++02Buxlm7Q1vjCdXy3WGsyz79WtzRhw7nWjaAtBfE3YPJ4teieNsC3WkboyvPSUOT191rpqevAelDmTfnyYTVD5BpBM1TSCGvQNr2f9FOdZn3C6HlhtZ8r3ouNrrQZQJyKM4G3fQDGlr3F17Op199q3n+xO1vJP3+afUOfLlki+nyOtQty5woSln7nqbLUcPjUe6m788zNkekhtq6AxGs3LQd2LMJ+Vt/Ua+PhfwfyTMWxkKJxB7W1WPfXkeiJqNzxoWSJztS/JK7UkRAtKQ6hATGld5U08DF/YyBIRgK46WvVuH2DCm5wpgt5pq78yYAeFX1fElod9LvC8b8BACo5hjfRIkg86DPfqXc8/V4SHg4ILso3LoAeP1clGimcnqyoalfnOPVViiygZzT5Hl4VsxbaX3IZwbFidety4z7gLMwJi/chI+mz9SVvJJ0b8mriVNdk22an/Q7TMm6Sw/aiD6kbGLiSRBi5muEeyzZOPDTs/nGngZ+2zuawSnHArW9DKySc0J+dvVOy3k8lZesJSZ2sjroVwctd876GuCVYQh9dKPhZle+ZRa+KbhDOxH34SFDn3ZNOtrHIQH2aQVWbEze4SlMycJ2esFdmFowWng9RAnjpn1Gk/FyiBoaydc+oRbWIcg27LI+eWqNlOex+CN7SpU2eX0wDu6DVZTsnZG4dP/rMuOeBHpJHYTU7Yto2HtXfBuXx3rEyGvfrsbWNeo2T5F4fRE/8Rmhg/B4i7tsKTMQShZ2IY91X/WNNWI8HcyQ1vWjjeCpvNc4FWCfkN9amzlEg9PGmrSJ7feucKQevMg5IV+6fYHlPFIbWXtqGh3sV6FtTbRC6pe8vjFiiPUbiIUPAHsbGrFo4+7MCXl3bsaALodh99A/4u7A9bjy8G5889dbDbJuNvn6i4nWK8KNWIsUI4y1NPlgiO+kbEPA/VEunTLXiCo154Q8D3htELHJ4JLoVJmb7OKov04P5r2N4Z7MK//4NyoGeDxoecrDmHr/2bh+RHfuJeiqRYn1YFLdd/3AoSZ8cHJT7c8W8V1zoOeTUeBzViyFbNoJzi5yTsjX1old7Smyu6WGOE4ieiqf6YtIeKp3FocaGSPR3Nq6pEDXRNIOlJsq64lXxuL2Bx+OlZxcD0+TJr/+4ydN5a/Fdk8bxeOibfKi2h5ptLtMmrxRO7KeZ9StUtxoWQ/O7ecrbfK6OPS7q3Sl+zHUz1T+Rhq10U6vJ+cq/3JDeVov0Uh2xvOrRQtTRd239Q684HkegILwSwgM13nuX03lr8UOr7KQz1YYGBhjmLJIYa/UTO/UaNA+91tr0I5qnK4CV3JOyOvFbFsz8pUvh7FYFGr9pR12QshXnrNN3owGFAxb109DU+5LrofK4iFjuGfIziOyoiYM+Hr5NtzwTvrOavtt+cZa1iltLOz4tGrzoxkLeXONzUh3G+E1FtQoTawkdJCW6/RF9zMCfw9K4xn29GzKnCgDB2JFcjUEC5J/l90W/3vdDasSyhVDfaBpbqYhEMKabcb34tWCMYYd+xpxqGdx2jlfULusjK88faLJYO2aD6IeTTMW8uYQKj60Qg1vzt6NhFN5OSjODRUAwMG7Ro1l4c7o0PuQ+O8ubcqaihVUZo+ED+GmXXU418d7FSZDu63f47/5Tyie2+c3vzipOcv4bdctdLoKAJqxkDeLyDaaulo38VdNnZhY03zR93RCHJtdrUIM7lnVuyznqzYZ6UEYfzy+d2LK+F+i2sZ+1HQ/3tqNQsoobFCO0VK9ow4HPDRVI/ZPRlUet49vWpGs5FxQf9UM1Bx6T8IRsSOxE/1PCc2/CXeYpoQLeSI6iYiWE9EqIhK/7FAnPq85bU/LE8Eq6TknbF4gpL04M/GaaE7ZysotFTnkifQNuf+vcbL2RZNvNV1egRfwehJehu1B3fnPFTCmbuDatH0XHvK9iTkrzIWFYAAmLtiU9DuVFl0OQfmJNoqGtgdYuvyL0CGZE8EtIl6wkCciL4CXAZwMoB+AC4nInFuLBo0GFwsBgM9r7tZFavKsbcqjSdDsRTQYUcuoMzGgY3n872Xhzpbyagikv/tOtF0hZQLz3jJdXnoYWvu6MjEGJsI9hakvhzrBPw1X+qbi6E3Kq2sz1oYxnOL5Kf5TX2wpse3y/RutxfxvK2y/1+x0oRwKYBVjbDVjrBHAeABn8C7kwznVvLNURaQPbVmLfNVzRneZcgK9VSzK92E92gGwHtSrB4kxX6jhSdWkbd7HN1Uc8wjIF5nQVba7e1hk1TWZHEEwxvBK/ovx3065UCaOGEsL+W3RqIVbPIlEt9COABJDOm6IHuNK2e6lJq4y610jrpV6NVz/RC2z4ZudvvwqywrjS8et2uenF9gTHydGupBPtMnb0amTy6jZaT3eOUWjlmpjbnuq716+3khyYSxkHOP869RmUkOWZyJbwxooPY2keyGi64hoDhHN2bbN3I4+HVoaX1BjNiw8p4gHinRuZW5hkFn4m2v05Vde3DTcZTZrwvpRn3hNwnYhn1KvoDFBopgjY/BlCliv0lYytaGjaz5I+h3W0UZEdDGfx36t2uhyDVE1FN3DNgBINLp2ApCkMjDGxjDGBjPGBldWVpoqpKTAzN4n5h5pzx0zTF2nj1TvmgSBkgXmGv1Q3OxVXmzvh80qZs0WvBClCf/JN0HxeOw9ae3kZAiTu1Cd7X8IX4UGmi62VYkoO7o6eazR9jKVEC3kfwbQi4i6EVE+gAsAZHB9cDflDXy2XDNKdoh4/d41MWHJI2KkEFSkqVdDSIm2RBzy0y0oW/a/pGO8vv1qG3jHhPzBtELxvFH0xP5SuqW5rA++Ch9kuly+/Uef2Mz3G9zbIhsXQzHGggBuATAVwFIAExhj6cvqrGKqpWeH2IxBIurLPayBPhorD4inLS2yX8PSy/LNtXj5q1VJx9LMNTbTMkXIFy54w3Ke1pqBsYtDOgpTS/FW6ARDZYmCtdPngqk2x7YlvwvP6mREuEGUMfYZY6w3Y6wHY0xpSZ11TAh5J8O3qpHW/hN+B0WM07lnqS/DQFlVXFi2a1nMuxLcOOuV7/HM1OSAcLq9qwrKMqcxQepiIk/NWh65qp7JtC7EaBNKze6uwHXpZZIHR/hfUKqNwdISyuXU3ycd8iYOulifGPOoyaU2vbjURS9unfUSjitN3Ak94OXpy3DRv5qCQwWEbOTsfEjVmp7cPWq5cOCOz1Cv8NBjW+7tYEobdkcb1cUfADeKiUe/uz55he++Rg77odro8sJSzF01LD2sMIFhysOXJh376KbDcPtIe4WjEgOGHQvy6psD9KpM9iaGwraDHBHyxiW21U1zrUwC6WG/r+/ElwV3x387F+NaP/pXAxO+LooMvQPFHQAADcwe32W9FAV34ySahUn59ycdj41A7mz5d1zfeHvSuTKK7k3aayRQbm2Rlxqpj3jNdr7BylLxB/lqF3r95FOdKQ7uUoHbR/ZWSZ0ZXkpdhQHzokdtwaXNHmW5IeTNvEGPGY+cCPtYAWaH+5q+Xg9ne78Tmj8gIkSD/vxG3vQPzDhnEdZFozGHXdgUX8p7EQM9q5OOxTT5v994Om675U7b65Ta0nksuNFa+5FRuTC6aYiu9PyH2VvbHal6bmlYv41c1QSjgE/FJq9mOpIbeWtgpkkESX11qR76dWxp6XoltF6yUsN4/+jplsorfXWIpeut0KqkAMf27xwPceCW1YGJKNUotoirZYs89OuQbHcPm118YYDUFrLXb13T3rFEPWZ8pjsyKpdS27iywIukmRg6zFjmCmxnZfjpoqU45PiLVNO0vGO2/gxTTDALw91R1TAuLdkv4W7q+xpITd44Zpb8hyy47jEQWpcWmb6eF+3LC52uQhJ6RwbhgqYPZLwjuHBRlNJev1reNY8HLxFZHUV27rPui318/RTVc5k0eaMrwI2ENXgmcL6hvJXwEjC8dwd4NBZDdSjXv1YjNZs1LLKf8Opw8r7CdwZuVNX61dyGZTx5DYS4F2YslH+ZRt8xuVAw6oFK2jb9HX2MPMMPi8Sr5UKZb/+H/yyPslnvs9BQPBy4DP0blAOL6aWvx1z0STX0OKDGPiwjerlvm8VUwd26XWdMuP5Q5KfY30/uv5+6aUdq8maw10+egVBRYrMWrVBdI/ZBO2BFFbrSVbVpcpukLGuCStp9jKN6t1U9x4t2CbHlAfX6LAp3wxuhk3DB8B5C62N4I28Dqvx1I7obrY4qvNpZapfzeQhDu7VKU9AqitU3slezyYvqztnVw1SwW9aVUj36ttcn0IygOTRWOOV1kZC/tfEWsDzzPu9uXLdgFKseW7y5dHhX3HlCH7GFWJx41bqax0RkU7viYwtJtRpoTVlTG5UPrJrwlxOvvLHYIYta8alGAobfscuEijn4RKN0A25ycj2yZ2s8dmZ/p6uRhtLOUGlE04RNxrkRSWqX+61saOS4wtv3epWdO+xWaLK/ZwHmVrxaFZD9zzZ96evBkwxfU0r1acccCKynidLcyKauZ2pfE70kFzR5N+5fKloPMGyu0fUpZNG8TVTIJLc33qQrXerz7DTkNMV0g7u1gkfVu0bNjCM1eXVMtGTLbd9C7wmobJSh9ZIP9qxKO6a62MIBGEgxtGqHK9/Erru24ZezlTefbhpM54CQ15FmXdhcpFXziH2uhr1r0vYxVq+fLq0/A3pzmBg+Qle61HmwY/rE5mGSS2pTUqAq5KUmbxNOCpVDurZWPP7Nsq2Yv26X4jkl3ORd069DGdqWKk9GVxTno12ZykrB2NA8x4X8BhbxFKlBCSZx8P+2Vht+BBS2X9RCj9zuGHVptCrj3wkeh4dKHo6VrJn2vRv0bQmo1ko/8aUGTyMNLxop5A2TbR4a+T5lTf7d2etw1iv6Y564yVxzaHflD1eM/AzLEnJBk9eSIx+GIpqi3eEpEhXPOvCP+HnE08YW5OlxrunSqkh3Wi18p7+A+6+NLoLKkNeQKn1zbGqT63O8B6anU/WTly6UxjFjOnFw0pLXfq0eF2nymahoqeKNFHcBzAEhr0HsI1bg9aBNibXV1mqMDx6NFSwSM6esIbI3T+I8SSP4xwfqz1YaSq9rMvWIOwDoC0usxQVDu2C/ltquzvU+YyvX1VrpdUcacfeUQt4w2eZkEvClR94DgP5Uja60GVN+/V1XPq7yk89Ul9J2ysejHdmFc5a6iZli9uWpa4N7WcQEsdbTUXUkZ5VT7n8fOzqNBADsKWwvpIxU/ptvLHp4qtxWHMFF20psUpfLPEZKMVc03o2z/Q+hxf3GFnupNfOK4uQPt6Yip3JKlAul+ShdLsLcUN9BTV6l7BfyXwEAjP14rq583CTjraL2TKrD7VDl2WJzbYyxMtwRnbzbNT9U2/M74CL/vSjoOgzn+D8UUo+ywrwmO3BUYOTlNXVxN5jEtCZTWUk70N6mdx2bpBVR6xPPvDTyvAyidy1ELN0mT3t0CCcrbUxFt95e68d+imeskRNCPtvI1NlOCHylKx83Lb4xW5NMo5G1rB2q4G4hr4fTRp2FxsI2OLR7a/y29QBcNqY93sp/SkBJkecZi7FDXneFcN64q071HP1hARBuipdvxVxTxwqQFGQipZ1dONT47kxhpj77l9oXY2YyPymYi1Ta/P2TFmFiX/7RbXNCyJuKXeOkfMwg2DqzTZrndWaTFbjK5GSRVJl0SeM92MjaYA1rj38VtsZJ/SN62qCqSrz15D3Ar72AD67mWof6vIiNucGgrdkuinctReLUQNLrT4n9U1poXjylTXDr+F5kGjUa+uREb0zpGjVNfnNN+loYHuSETT77pB2f+voS3GsOafgnvgvp23vSTWRadO4GE4NeUu/hu/AArGHtcUTPNhjRW8GuPOAc7FTYGckKNPRajA5cAxp6Ddd8efF43n+Sfmt95NuURLyB1NYU8SZzxE39bZEU/orRUNRBf6U4kBuavCk54OD3jdNHyedtmsC7+NhDUD/TudDDZj2GMsYrN5WrM6Rq8uOuHYaurYvjft9K+DhvDH5M3/Y48J4n0brEvRukJ6I1Cm96nvw+9DtQBjVn3/1op+a12kKeNH8mEshXHmWJUmdyQ5M3g4MKYn6ptk+5XhJv4Y4T+qBjhfMx7o3iSZ4nTCN79Ph0DuvRRlPAA9rx6c2iLuBd+DQ1qmTF2yR9cR1lKg4FZGG/XBWbfEFeuidVUCWQ3x9975svX4PcEPJZ5C8OACFfEfDAdsv5ZJ2VSoHE4frAhjFp54sVYvY4yfyKEzXOGhdKPnJfEC47CYTUn1k42q9rqUko6t99Sy0+jHnM2OTbt0z/yBdEXWhnFySvsj3P943ZqmmSXdKRI07KRwYCOHg9pHnXZNnHDgAKW0Q6cIsO/fDhnaPSzsf8ywFgaVjM5thG6H7JP1TPmREgIjT5TFzceI/tZapRSTWq54KlnfFo4FLcnddUX614/nrQsruHMnxADM0PRfum0o5UHq+YdRJqZJ9UUCCbJucAfh8YSm1Ajqr25sr2lO0HXDYJrS9/Bz0q0ychN1WdFf87pBLYzU5atlZZ1AVzsVY0d5riTKx6L952sW1lZuI4z3zVc8UFPrweOhn9+/aLH9OryafJBB19Yw3LtIDMgE0+xmF/0LjCnhmnnBDypqIa8K+Gbni92rRJqyzU5AEA3Y8GCsuUz7XsBCDSuZcz5zV5LczYkH8sTQ1sJY6Y4HNTSGStqrRskYcfRh+LR04/QFd6vfnaQXlRdAXsgHMUzsbmB+z5wGepVEgly6Q8J2K3sDe27CNbhbwGngR/49/C9rqeGcWMYAl1GMy9HtlEpuijHcpbIC8hpHZNz7M0Uifmy78vGHm/eVp+n9F+atHypBtLT4KIniGiZUT0CxF9RETlCefuIaJVRLSciLRmqyxjZuVngJxzMeNtrokN++yObieCdaeNT/qdeE9uN8uZs8k3b4y+U+/pL+jMVxmt0koyTPJrfTjS7qNYOd7O1u5nJ7gbZ4e5ZhqA/oyxAwGsAHAPABBRPwAXADgAwEkAXiEi5w2qCfCwh23pqrwrjF14UlfVOSrk+QjgUGFykK/ED3guxJxPxa4hezJOGzOaMCrkyWPWYSFzOY2Zlg1pmuQjJ6s9XYGHdwN56WtWwozQ9rLXE7JLfg+/5h+UsY5msCQVGGNfMMZizqU/AegU/fsMAOMZY37G2BoAqwAMtVKWFmY2z+AhLuoqB5q8ko+wSsvFQSEvas43Ucgf2FklXLFbMCE77dwTwI0jIcNmFZ23oHavWsrdXqa9zkRrLsNf3gMzQgfh+eLbNeoUI2qu0SyNHzylwlUAPo/+3RHA+oRzG6LHBGH8ceXnWV/s21AWiSG9KN+osOekSaVJVvd1Yqs0jVYI+7V092KvbHGhtLwbB0fE1cR4X8ikqGh9JNu2LMFVgbvRqb/6rl/x68neideMko6IvgQUI2DexxibFE1zH4AggHdjlymkV3yfRHQdgOsAoEsX45HhzHLgiDOADyyuMIu+rCA5Ex0iLe5HDqyOKi5IGY4n3BOvzVZEYca7xs47stcSrA/j5hq9LpTGsSTkywox+77j0KZYfa4v9XrXuFAyxkYyxvor/BcT8JcDOBXAxayplW8AkOjv1gmAYmhFxtgYxthgxtjgykqTmwOY6CktBpxuriyFgslw5+Zrrok3Fo1pj32F7YAL/sulXJG07ZocajXRXONGc0MiZlwTa4vc7RbqNvS2ADU/eS3BmkkYZnq9bUsLFRc/pV6fVRt5E9FJAP4M4HTGWGKg6MkALiCiAiLqBqAXgNlWynIlTmuWKW5aWt413/V/HNj/FGFV0dNwz/Y/hPsDV2onymuBsK9plWvSI84B76FU+o4427ay3KTBxzA6ma7Xky71XkuiG4QU51sZdZvr77FIo3WITsbGXChTailKs7faa14CUApgGhEtIKJ/AQBjbDGACQCWAJgC4GbGWMhiWaqY3Tzj2cC5eDd4nPXyHeo+5Is0mtdxZvSAgxOvOtLccsUlaHX0TRnTea74tCnfLPKuMdMKMgUw48FFjffifP8D8d+itpkzh0EhrzN9qtIRU7B9FrqI2ad2T+BaAMBSVhU5EAt5YJPcsGRMZoz11Dj3BABjG0CaxaSQP+76Z1DbEATGqd5GpoKj/3cmyJTHm4eqhnEoK/ThVsD5kUUGjunTFsf0aZs5YadDsCzcGft71iM/Ic7Hok21GCawflZhJjvt2ODJuNoX8VnYhxYIDbkWZT+/yK1eA0ecgTyvB+z7l7jlyQtRpgszYQ0ypTBb1zZt2gB7EF/UlVXmmmzn4C4Vyps56MVhoRpT3GN2wK2l6puGuEt7y0xM4zqwU3n8WN8O5YppXYPJZ3zoTWPwj4MmA4gM6ev6X8KzVvjzSfvjjuN7x3+b/RiJQH9UyQh6u1zaHRa0BAacC1w0wVB5iZTTPlPXXTC8R6RO+ZFgfHWt+2M3K8J7xReZrosRcmPTEKe+VfFNk1XODzgP+DW9UfH6kpcW+HDrsT1x6oGR5f67SjRGJIKFfL5fe8MFo/g8AMJNnZqBkjZJcSNmH3G/DmWoGNwFWMC1Oqq46Xtv2LtGp0kyLV+PBzj7tQyZQ8jExb52Q/DXwAVY0e4sDAEQyi/DQP9rGFJgz7qPHNHkndGoKZO/69mvKl/HqSUREe48oQ/67FcaO6CeWHDPrt/4q5iMo516IeuBkNfdux1xe8LCR4jukfKG/dJ0PxoXmS6J8M/Q6aj1OLPvbm4IeUHvM+TTXnxjagPxyIVCGNFb3d4tulsHQ2LmJcjjxRn+R3Fl491oNcye4a1ZXCRWFHG7CypPpocHGb5G1NMpKYgYTDpFJ9lbFUc8ffq1V4m8ypkcEfJiXs+OS6drnq/Y/wh8GhqKVUMfM5QvM2iH1Euhxipe0XZYsx5OqkRHHgyEhawnalGEPh3cGdagiPwArA2W7BC/8cVQLrLXGDfX6Ev/aOgKE7URQ/+OLfHPiwfh8bP6AwB6ti3FRzcdhvtG9ctwJR9ywiYvCpbBPNC2ogwnPDwVPg8B39tUKS00wpsK79iCPrTcPx4CGOZZBgDoXPeL6Txib4ealb4N9PWsM5Re7+j5iN5Ki/T58W1oAEZ49ZsoTx6QvCHJwV3SFRZRrtg5IuSd6xaJsa71IkpupYU5SCAs2MuT9y3FV/NG7+muE/twLoE/ZOUh2/kxc5Emf5DnN0Pp9T6mly82bq4xwpWBu5EXCGIZxzxFfd6luUYrWxMPfe3wRwXURB/kUfc+YUyslOf/Blg0X0L1U6Nw8zFm1zKIZz4ioRiKeqgHp8oExTetLoEopcWNO0MZRe/IrjDPuCeWkaceghcN4OsIIDV5DdReTj3LRwtqtJCx8c5W2/YQ8+VZRqu+WWaTj2csJluuFJYBDUCLfPMunuGiSjwYuBy/FB+Gf3OsWiIVxXnAPqB9S3d7KWkh0ny3l4p1d5OZdx+DjTXam4y4hZzW5K0OfwL+usyJUujTrtRSmVbQitCXbdpbRYuI/lFRnO9wTTLDQ0OuLC3AvP3OxR3nWA+zoUZedJ2Bm+Y5jC6GEgmNelZ32s6tijC8e2u+5bs0do0rUGu0Vh9ZyIRbYF6FnqiCosxLzk288hYcZS0ibmZF+WZ3ArKR+IIt8+R5Pfjk1iMjK7BFyb1z3wD6jAJK22dMahc7YEwpEvlJaNNGR8gNgfRrFLPWJCeEvChIazNeNYpaZU4jCK36hgULeX+Ic/6n/wPoNBQo78o3XwHEJ4l5PWNRmnbnocCF4wCNuRu7qYOxAG0iByFOBRoUTU4IeVHDTy1vFTeiuaGCYCG/rbaBb4ZdDwOumQb4ssdcIzGOV0Coipmh/qauy7b+rpecEPJqnczw/pEpUFhYdGRBaHnXiC05N7uHPjbnRUx0/kJ3LtZyM8xnVJPP3NLK25j0kbdzw10FjG8jqo+cEPJqlJL+2e9NTMHMkgVaZCKau9KIdqHMUS1IDx9VXIWLG+/B7tZ8fLNNh8swwLywO1xSdw/KvL+AUTwmt66QmryL0bvvoxoDGl7Dsf6/KeTr3sezuyB98kxL0IqOeJ+b3UMfpw3qiu/DA5oCxWUBN3gfwYENY5yuBna35q+9ktn9iTLEqhJN/8aFQvLNCT95q9RC+eXqDWsKAMf7n8YeVoRZvCqVgTUVh+OglGNa9T1AcDCkXJ200sPpAzvg9IEdOOYo/pP5yR0jsammARgrvChNGkq6cM+TTI5avSV8XSLdQk4Ieavx5H8YfWxkhdwzycc9BjT5i089AV1bF1uqhyEUNu3WGnnEIt+JI7uFfOLuTM2BtqWFaFta6HQ1hHjLmBXyrbJgTYYZckLIW9V8Oqjss2nEznzF4d30F8ijZSvkoaXJZ9tiKLsJ5EpXkICQbQ4TYnGv0dkFGDHXGKFlC+sCRfED5OBG3l6bjPL1hWIWrLhpk/Acnf9TpGM5fzu4WU0+V8kJIW914lU9YzH59hx2quU8mJJA16ivaEXeZ5P72dK+twnJt0vrEiH5mqP5SHkRk9UeFuSeZzaTE0JeWJgAAd4100MHAwXWBYriKIOzJn91451c8+MB84qxm7Yrc94+LeGD1OSTyQkhL0rvscNf2TSKWjvfAGUH9u6lvzrGszeJoDC8Dpq6JJyRQj6J3GjZosIauNhPXlFr52yu8RowtPvselTCTHPuedcuVi2ygj11fqer4Cpyw6VAWL832fFH/Q1orbyikJs/uUFzjZlSV2zZq786ts0WChLGzWm2M9eRmnwSuSHkhek+JgXKkGtUT4UqupusSzKGg2KZkPL+oAs7iyBhzDTi/kiyixJkx2YedsFFLSKiPxERI6I2CcfuIaJVRLSciE7kUY4awvYdEGAamFBxLZ+MDJtrjEv5POj3UrBtxavCIjAu2bpIkzfzriRN9PJsdLoKrsKykCeizgCOB7Au4Vg/ABcAOADASQBeIRLUO6E+QVrntbaUn0fHv6zxzzjJ/1T896P/d3D8780wv4zabHjbeS2P1512UHiR7rR2yUhR5bhr4lUKeQk/eLTs5wHcjeSWeQaA8YwxP2NsDYBVAIZyKEsF5Z6/pJW1rdR4yJNvwwOxjDXF52jfsml1rSXt16R3jajY5672RNKBmzR5uTxZm2khJ/dRzj4sCXkiOh3ARsZYavi0jgDWJ/zeED0mBLUO2uhRXk23j+nbyNjNUSjbDFDQyHUJKv3CzJDcs0tGCpJ/btLkpYzXZjnr5HQVsoqME69E9CUApSj89wG4F8AJSpcpHFNsukR0HYDrAKBLF74R6dQ6rt7NRHhop9PvPAqlBT7gOctZJdH5gMMUjvK1yXeqKAI260trm4znKIz3sQIUU9TdzkVCXnxg6OxmGyt3ugpZRUYhzxgbqXSciAYA6AZgYVST7gRgHhENRURzT9zRuhOATSr5jwEwBgAGDx5sSocxOtTWa7LgES6hR6X66lZruSvFrtFjrtFPUb5+56sGrz1hAcK9TgB+FJGze8w1LCxVeS3k0zGGafWFMfYrY6wtY6yKMVaFiGAfxBjbDGAygAuIqICIugHoBWA2lxobqqTy4aBuTd5N2l0KJm3IM1ufo78IA91pr8+eDcyHdOcXoCzxCTI32eSlGNPkuL4mt/drpgjxk2eMLSaiCQCWAAgCuJkxs9u16EGlg6ocDukV8jYF3QqD4DHcsQ1OvEbPeVv3AFbrLMHA7ZcX65vncBOJHzFXTbxKNMnzyTUNRuCmqkY1+u0Jv59gjPVgjPVhjAndjcHoBKn+Db7FdvyYkAl4TATHUgw1rFHfqL3m0kO76i6iPr9N5kRR8u2KNSwI5qJRm9TjM2FPW9vd9yJbyhGNe1q2jeieeHV493ZtjGryEYzc0vpWh1qqTTaxuv0pTlchDgvLiVdNbBp1tTz/n7aUI5qcEPJqXjBqGpHeDSKyziavueLVRMcw0JncZdPWR+IzCeSV4bOQwKUcEm6I2swnV8mNp2VQwISYOzR5a6EA9Gnyu1kRxgZPxp4OR0QOGCiyuED/lE02LoaaUHpZ/O/zBmd23w0Li5+RgnSU1yYLFQonyQ0hrwKp9BXdE69uFlw6NfltrByPBS8FeY3PsR/RS78nS5s9+kMguIVunTrE/25R4ENlSfZNHjdLMgj5AJMTs4nkhJA36hmh3ybv5sejT5P3eAh/OK4XDu5cHj2iX0s08lxb7VmmO61rcc03XWry2mi/qDySG3knkhuhhlWEkdroWq8mb1evN7kCTFcyr4dwx/G9zZSQ87ZPUhvqqWCX6JWLobRxt/LlPnLiaalOvKr0Fd0rXnna/v4wH7jpp6RDY7znAwDqPCY2MzbqQhmlrNx85EtNcsBOalu45Iy4pR5uJfvbmp3khJBXe+kbN21QPB7SuUEEV7nVqjvQtm/SoXVV56OqYRxCxGtAlbnChrQgV7uQWsetd8fkzkaayIVrxsgRc43y4Qu8Xyke1x3WQHBjev78g3Dz1r3wvsUpQ+71NRCxMgs9QtI1d+37FRWmWWIQKeQNkRuavMGX/lO4n6CKGKNFvhcDOrXkKDz4Nv69/uQJrBpWzDV/pyEAI/1P467AddHf2h8q28w52fe9tJncEFt20Syf1qvBUU5XIQn9YRYywFnDqU0R8jPCB6ukBLJTMjFs8HXBe6Gjo7/dosln47O0D2muMUZumGu0+gR5gZTYaEO664/JYgdBXptIKzR+K92hRcPWpN+5JnoIwBe3H4Wlm/foTG/PE2BZaPqyFSnkDZETmnxYpVPUecuAPy5OO17VtlxwjYyxl7XInMgkeieZlWD5yV4/hSUVqmn3FnVWPedWgr4idGldhBMPiISurfWp3x8AbEO5DbUCwl7ji7LqWb6AmhjDrkVIYZ+4/pKL5ISQZwAeCVyKxwKXAAD+EzwR68KVeGngRKCsPcYFjwUA3Be4CgvD3XHLKcp7RJ7hfxSzwvsDQNxO+0HoCNSzfDwbODee7orGu7jW/+GWjyge/zLUZB5J3NdyXrinal5fhwbG//4h1A8PFt2blubwhr/H//42NCDp3KJwFV4MngkACLfsjEMaIkGafgr3xUetrsEdjTfE81h+4N3x674f8qJqnXhze+NNacfqWT7uD1yZdvyJwEV4MnBh/LefNQ1eq9ufnJT2/dbX4+nA+fgudEBaPrPC++Ns/8MWaq2fYMuu8fZ3mv9xvB0cCT/LS0ozyv9k0u+D/GNMlXVidJP56nA7jAseg+mhg1XDN/wQisxlPRa4BBc23oddLLJRzMeh4bg/cCUO8o/Bb+H2eDl4Ot4ORvYaej14EtaHK+FneVgXrgQATAwdhqcD5+PRwKWm6ryl/TGYH+6JRubFjY23xY/PCUfWgwxsMPcsAOBI//NYFK7ChY33xY/F7vNM/6Om81Uj9vwBYFjDS9zzB3LEXMMYw39CkQ47NhSJJvgILse13sgX/97gNbg3eA0A4N3QSCxTGe4tZD1xfuOD8d/PALgzcBPuRESovBQ6S0j9A3llqGoYp3JS+e9qlbyuCPw5KV1PT/qOTRtRmVxeIC0JnguehykE7EDLeNojPS3wYXgEPmwYAQBoc8KFqJp9EADgaW5uoJmZGD4CExuOUDz3Tkhh71sAY0KnpR17LGWxV9DTAq+EzsAroTMUn4ldMADvhY6OzxX8GuyOB4JXpaVTbTMGWM66GMsn4bkcrPBhOa7xb/G/Y3V+FJelpYvxoOoZDciDsxqbBG5q/a84rMpMrgCAG88ciVM/aqeaN28Sn/9lBsKAGyEnNHleCwSL8p2JefHEmf2F5W3Fepnv1W4erbMg1sudx/fG5FsOVz5p0LZrlym4VZHzphezPHfewIxpLrcozPp3bGnpei0uGsZ3n2kjtCszsa+EDnJCyBvdpNqnsshnzv2K29kKp2OFu2yMVa2L8LdzB6J7yv60Vx5eFf/7hqN6JF/kgrmwR89IN7MUFfjg0Smd3TKfV1GcjwUPKo9I1DjpAOe2xGtbGvnYP3Zmf/zfoE747Unt2PyPnNEff7/gIDx/fuYPghI9KksU33WMAzvx/whUtS7Cl3eM4J5vIucPETOvlRNCvrQgT/G4mpOCT0VDVdq4+oFT+2Fw1wr8fN9IVD8lxvWyKN+H/h3LVM+vTuk0+++nHgZhSJX25KEWFw7tjE9uPQKTbz0CZx/SKelccb4Xx+7fNIy9+8Q+AIDj+0WOOS0fzx7UCQd0SO7c5w3uhEuHd0XX1kW68kj7cDlIeVSb71iuTwGw+oH60wm9MeX2IxXb+JWHV6H6qVE4qndl2rkBHVvib1Ht/bj9I1FLvR7CPy8eFE8z9vLBAIDOrVrg6bMPBACccVBHnHVwJ5hlZN9Iu5t8y+G4OEX7FvGx/vquY9CzrYnwIwZoI2hknBM2+Rb5Xqx84mQ0BsM44KGpqumGd2+Fn1bvNJT31Ud0w9VHdIv//vCmw7DPHzRdVzXeumoYBj02LWO6OfePRLHCxyjGf68djs17GrDs91pc89Yc3D5Sf3CyHpUlaUPh1sX52LGvEWOvGJJ03BMdDXVrE1kgVVqo/KEVyaAu5Zi3rgYA0LNtSZIGlyis8n0erPlL5EPZ7Z7P4sdTZcHAzuW4+6Q+eHrKcmF1NsIbVw5Bv/ZlGP6X6RlNkg+dpq7ZajH65P3xwdwNuPmYnnH/8/KiPNTUNRnfY3n/65JDsLGmHiOf+wYA0L2yGPeN6ovh3VunfRxOHtAeVa2L0LtdKY7r2w7fjz5W9wdLDx3KW8TLnDh/E4CIBv/Lht3o215dYXIbk285HKe/9L3QMnJCyANAnteDPK8H5w/ujP/NWQ8AOKJXsj/8u9cMR0inAb/Ap6ztD+piXlPWolVxPjpVtMCGXfV4+aJBuHncvKTzn/3hSOyqa8z4tfd5PehUUYROFUWGRx5XHd4t7VhMK+oRNd0seuTEpPN3HN8b3dsU48QD2qVeKoxzDumE9+duwJBureJCvqTQhzyNOQSlBTRKGt8NI3rg/w7uhEAojLZlBehz/5Sm9JZrboyj+0Q049n3jcTehiB21TXirFd+ABDRoH/duDuedr+W5uy5NxzVI20E89FNh+OYZ7/G2MsH47i+Te+1Rb4XPds2mfBm3Hm0Zt5f33VM/G+eAj6VPvtF6nTbcb0wrHtrlBjY7EaJ9244FHsb+CtyShzYqRyA2DBROSPkYzx+Vn9cO6IbOrcqQkHKru5eD8Gb4Wm+f8OhWLezDiMUhqai+fpPR4Mh8sFasL4bPpq/Edv3NgIA+nXgp50c2asNZq7cnnbco6OlpXagwjwvLhhq72TV02cfiCfO6g8vEU47sANmrtyOC3XaM4vyvahrjCyOU4pe6vFQksAc2LkcC9fXAABeumhQWno7aFNSgDYlBahCMb744wjUNgSxu74RV70xB51btcA/L1Z2CTZLtzbFwkyTIjhvcGcc0KEltwnZIVWtuOSjl3kPHA+fV5yUzzkhn+f1WLKdDa5qhcE2v+QYiXMF943qh7tO3B819Y26hK8RXr9iCPzBMPonmLZiZpdUzh7UCf/+drVl7YgnHg+hwBP5gPfvaKxzL3n0JMxcuQ2Xjp2Nod0yj8rOHtQRC9fX4Of7RqKy1Hlvot7tIm2bMYb/XDEER/ZqozrHJIqLh3XBu7PW2VqmFkQk1ONGzVGDF62KxXpTkZuWUA8ePJjNmTOHe75Voz8FgKzSTuxg575GBEJhTdetcJjBHwyjhUPupUaZvHAT2pUWYFh3PnHzGWMIM2QcAUpykznVO9G+vIVQcxMPiGguY2yw0jn3qGcCGXftMGzZ0+B0NVyHHg3C46GsEfAAcPrADpkTGYCIIHAkLXE5To3qedIshPxhPdwVkEwikUjsIif85CUSiUSijBTyEolEksNYFvJEdCsRLSeixUT0dMLxe4hoVfTciVp5SCQSiUQMlmzyRHQMgDMAHMgY8xNR2+jxfgAuAHAAgA4AviSi3oyl7N4hkUgkEqFY1eRvBPAUY8wPAIyx2FZCZwAYzxjzM8bWAFgFYKjFsiQSiURiEKtCvjeAI4loFhF9Q0SxACcdAaxPSLcheiwNIrqOiOYQ0Zxt27ZZrI5EIpFIEsloriGiLwEoxTG9L3p9BYDhAIYAmEBE3aEc5kNx1RVjbAyAMUBkMZS+akskEolEDxmFPGNMNcg6Ed0I4EMWWTY7m4jCANogorknBhPpBGCTxbpKJBKJxCCWwhoQ0Q0AOjDGHiSi3gCmA+gCoB+AcYjY4TtEj/fKNPFKRNsArDVdocgHJj3yVvNGPpN05DNJRz6TdLLpmXRljClGVbS64vV1AK8T0SIAjQAuj2r1i4loAoAlAIIAbtbjWaNWSb0Q0Ry1+A3NFflM0pHPJB35TNLJlWdiScgzxhoBXKJy7gkAT1jJXyKRSCTWkCteJRKJJIfJNSE/xukKuBD5TNKRzyQd+UzSyYln4qp48hKJRCLhS65p8hKJRCJJICeEPBGdFA2EtoqIRjtdHzdARK8T0dao51Ozh4g6E9FXRLQ0GkzvNqfr5AaIqJCIZhPRwuhzecTpOrkBIvIS0Xwi+sTpulgl64U8EXkBvAzgZET88y+MBkhr7rwB4CSnK+EiggDuZIz1RWSF9s2ynQAA/ACOZYwNBHAQgJOIaLizVXIFtwFY6nQleJD1Qh6RBVerGGOroy6d4xEJkNasYYx9C2Cn0/VwC4yx3xlj86J/1yLSgRXjKTUnWIS90Z950f+a9UQdEXUCMArAa07XhQe5IOR1B0OTSACAiKoAHAxglsNVcQVR08QCAFsBTGOMNffn8gKAuwGEHa4HF3JByOsOhiaREFEJgA8A3M4Y2+N0fdwAYyzEGDsIkRhTQ4mov8NVcgwiOhXAVsbYXKfrwotcEPIyGJpEF0SUh4iAf5cx9qHT9XEbjLEaAF+jec/lHA7gdCKqRsT0eywRveNslayRC0L+ZwC9iKgbEeUjsiPVZIfrJHEZREQAxgJYyhh7zun6uAUiqiSi8ujfLQCMBLDM0Uo5CGPsHsZYJ8ZYFSKyZAZjTDF0S7aQ9UKeMRYEcAuAqYhMpk1gjC12tlbOQ0T/BfAjgD5EtIGIrna6Tg5zOIBLEdHMFkT/O8XpSrmA9gC+IqJfEFGYpjHGst5tUNKEXPEqkUgkOUzWa/ISiUQiUUcKeYlEIslhpJCXSCSSHEYKeYlEIslhpJCXSCQShzAaSJCIziOiJdFgcuN0XSO9ayQSicQZiGgEgL0A3mKMaa40JqJeACYgElBuFxG1ZYxtzVSG1OQlEonEIZQCCRJRDyKaQkRziWgmEe0fPXUtgJcZY7ui12YU8IAU8hKJROI2xgC4lTF2CIA/AXglerw3gN5E9D0R/UREusJP+ARVUiKRSCQGiQbQOwzAe5FIHACAgui/PgC9AByNSIyumUTUPxpzSBUp5CUSicQ9eADURKOCprIBwE+MsQCANUS0HBGh/3OmDCUSiUTiAqLhr9cQ0blAJLAeEQ2Mnp4I4Jjo8TaImG9WZ8pTCnmJRCJxCJVAghcDuJqIFgJYjKad7qYC2EFESwB8BeAuxtiOjGVIF0qJRCLJXaQmL5FIJDmMFPISiUSSw0ghL5FIJDmMFPISiUSSw0ghL5FIJDmMFPISiUSSw0ghL5FIJDmMFPISiUSSw/w/Yx/SKxr/TREAAAAASUVORK5CYII=\n",
-      "text/plain": [
-       "<Figure size 432x288 with 1 Axes>"
-      ]
-     },
-     "metadata": {
-      "needs_background": "light"
-     },
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "plt.plot(cv)\n",
     "plt.plot(la)"
@@ -191,44 +169,21 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": null,
    "id": "b059cb98",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Making Boltzmann factors\n",
-      "\n",
-      "All ready!\n",
-      "\n",
-      "Start of the self-consistent iteration.\n",
-      "========================================================================\n",
-      "Iter    1:\tConv=  0.0367285320\tConv_errvec=    513.943248\n",
-      "Iter  100:\tConv=  0.0011664181\tConv_errvec=     30.442355\n",
-      "Iter  200:\tConv=  0.0007158388\tConv_errvec=     16.592697\n",
-      "Iter  300:\tConv=  0.0004416551\tConv_errvec=      9.952015\n",
-      "Iter  400:\tConv=  0.0002725440\tConv_errvec=      6.129032\n",
-      "Iter  500:\tConv=  0.0001681903\tConv_errvec=      3.789181\n",
-      "Iter  600:\tConv=  0.0001037929\tConv_errvec=      2.339308\n",
-      "========================================================================\n",
-      "Converged after 608 iterations!\n",
-      "Max error vector: 2.250724530684238\n",
-      "========================================================================\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "# run MBAR and compute free energy profile and probability density from statistical weights\n",
     "traj_list, indices, meta_f = mbar.get_windows(grid, cv, la, ext_sigma, equil_temp=300.0)\n",
     "pot_list, _, _ = mbar.get_windows(grid, deltaU, la, ext_sigma, equil_temp=300.0)\n",
-    "weights = mbar.run_mbar(traj_list, meta_f, conv=1.0e-4, conv_errvec=None, outfreq=100, equil_temp=300.0, dU_list=pot_list)"
+    "exp_U, frames_per_traj = mbar.build_boltzmann(traj_list, meta_f, equil_temp=300.0, dU_list=pot_list)\n",
+    "weights = mbar.run_mbar(exp_U, frames_per_traj, max_iter=10000, conv=1.0e-4, conv_errvec=1.0, outfreq=100, device='cpu')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": null,
    "id": "977b5279",
    "metadata": {},
    "outputs": [],
@@ -239,7 +194,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": null,
    "id": "53b9cee4",
    "metadata": {},
    "outputs": [],
@@ -251,7 +206,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": null,
    "id": "1ebcfbed",
    "metadata": {},
    "outputs": [],
@@ -264,27 +219,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": null,
    "id": "077d9fb8",
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAGoCAYAAABL+58oAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAB5NElEQVR4nO3dd3hUZfbA8e9JCL13BUIvIioIgoqCimBvYG9g3+KqP1111XV1XdvqurprWTsW3BV7RxQUERFBFBGR3nvvkBByfn/cmcy9U5KZybTMnM/zzOO9733vvSdByMlbRVUxxhhjjMkmeekOwBhjjDEm0SzBMcYYY0zWsQTHGGOMMVnHEhxjjDHGZB1LcIwxxhiTdSzBMcYYY0zWsQQny4nIiyKyTkRmRVG3rYiMF5GZIjJBRFqnIkZjjDEm0SzByX4vASdGWfcfwCuqejBwD/BAsoIyxhhjkskSnCynqhOBTe4yEekoIp+KyHQR+VpEuvkudQfG+46/BM5IYajGGGNMwliCk5ueBf6gqr2BPwJP+cp/Aob5js8C6olIkzTEZ4wxxlRKtXQHYFJLROoCRwJvioi/uIbvv38EnhCREcBEYCVQkuoYjTHGmMqyBCf35AFbVLVn8AVVXQUMhbJEaJiqbk1teMYYY0zlWRdVjlHVbcBiETkHQByH+I6bioj//4nbgBfTFKYxxhhTKZbgZDkR+R/wLdBVRFaIyBXARcAVIvIT8AuBwcTHAHNFZB7QArgvDSEbY4wxlSaqmu4YjDHGGGMSysbgZCERsazVGGNMTlFVcZ9bF5Uxxhhjso4lOMYYY4zJOtZFleVsjJUxxphs5VrPLYS14BhjjDEm61iCY4wxxpisYwmOMcYYY7KOJTjGGGOMyTo5n+CISAMROVdE/iMi34nIRhHZKyKbReQnEXlKRA5LwnvriMiNIvKNiKwTkT0islRERovICYl+nzHGGJNLcnolYxG5BbiHwG7a5RkFXKOquxLw3l7AW0CHcqr9F7hMVYvjeH7ZH2ou//kaY4zJbu5ZVMEL/eX6NPEuBJKbRcA4YAawAWgEDAKGAfnAxUBzETlJVUvjfaGItAXG4Oz1BDAVJ3naABwEXA00AS4ESoFL4n2XMcYYk6tyvQXnOWA/4GFgoob5ZojI0cAnQF1f0eWqOrIS73wXONN3+iJwlTthEpFC4Gug0Fd0qqp+HOM7rAXHGGNM1iuvBSfXE5xGqro5inrXAo/7Tieq6sA433cITgsRwDKgq6ruCVPvZMCf1HyvqjGNAbIExxhjTC4oL8HJ6UHG0SQ3Pm+6jg+qxCvPcx0/Gy658RkDLPAd9xGRjpV4pzEmFbavhe1r0h2FMcYnpxOcGGx3HdeqxHOGuI4/jVTJ11U21lVks6qMyVSqMPFheLQ7PNINvnoo3REZY7AEJ1o9XMdL43mAiOQBB/hOS4CfKrjl+wjvj/W9YT/GmARQhS/vhy/uhdISQOHL+2DW2+mOzJicZwlOdK52Hcc04NelNVDbd7xSVUsqqO9OpLrE+U5jTDJNeBAmhmmx+eA6WD8v9fEYY8pYglMBETkSuMx3ugd4LM5HNXQdb4ii/sYI9xpjMsGEB+GrB8NfK94Bb1wKxTtTG5MxpowlOOUQkZbAGwS+T3eq6vI4H1fXdRxpcLHbbtdxvTjfaYxJhq8ehgkPeMvyq3vP1/8KH93odGMZY1IuZQmOiDQXkQEiMiBV76wMEakDvA+08hV9DDySoMen7F88VQ37McbEackk+PJeb1n1ejDiEzjoXG/5zNfhh5dTF5sxpkwqVzI+CRiJszpvRq+gLCI1gQ+Avr6ib4Dzwi0EGIMdruNoZmK562yPWMsYk1pT/uM9r14XLn4b2hwGzQ+ANTNh/ZzA9Y9vgnW/woCboU7T1MZqTA5LRxdVRk/hEZHqwDvAcb6iqcDJqlrZzvQtruMmUdR319kSqZIxJoW2rYK5Y7xl57wMhf2c4xp14dxXoKBO4HppCXz3NPyrp9O1ZeNyjEkJG4PjIiIFOIv6neQr+hE4UVW3JeDxKwD/Rp2tRaSiVqy2rmObjmFMJvjhVdB9gfPm3aHTIG+dZl3h9H+H3lu83ena+ncvmPE/G5tjTJJZguPjSzj+B5zuK/oZGBzDasfl8u039avvtBpwSAW39HEdz0pEDMaYSthXEjqepvdlEG5dqYPOhvNeg4aFodd2rIX3fgMvnQrr5oReN8YkRIVjYUTkiwS9q2WCnpNwIpIPvIqzczjAbOB4Vd0Y+a64jAV6+45PAKZHiEfwrl48Nlw9Y0wKLfgctq0MnBfUhkPOi1z/gFOh82CY9oKz0vHuTd7rSyfB0/3hyD/AgFugeu3wzzHGxKXCzTZFpJTEzfoRnJ0I8hP0vErzrTA8ErjUVzQPGKiqCd9UxjbbNKYKe+0cmP9Z4LzXxXDGk9Hdu2crTHoMpjwFJWFWiWjUDoa9AK37hF4zxkSUqM02JQGfjOJrKXmGQHKzADg2nuRGRF4SEfV97g5XR1V/At7znRYCT/gSLPdzCgH3NI2wzzLGpNCWZTD/c29Zn8ujv79mAzj+Lvjdt9Dp+NDrm5fAiyfA1/+E0tJKhWqMcUQzXXs1TvfSh6p6RrwvEpHhOC0lmeQ+4Erf8V7g30DfKPZq+kxVd1VUKYIbgCOAFsAVQA8ReRVn5eKDgGsIzKB6TVXj3RrCGJMo01/G05Dd8mDY/9DYn9O4A1z0Fsx+Hz79E2xfHbhWWgLj/wqLJsBZz0D9/SobtTE5LZoEZypwBhBTN0kVcaTruAAnwYlGe2BJPC9U1aUichLwFtAB6Of7BPsvEMOviMaYpNi3F3581VvW5/Lwg4ujIQIHnunMvhp3N0x73nt98VfO2JwL37AuK2MqIZouqmm+/7YQkdbJDCZXqOqPwMHATcC3OHtTFQHL8U1TV9WLVLU4fVEaYwCY+4kz88mvel1nllRl1agHpzwC5/8PajX2Xtu10dnLau/u8PcaYyoUzSDj44HPcNpnz1bVd+N6UaCLKqMGGWcjG2RsTIKUlsLIE2H5d4GyPpfDqY8m9j3bVsE7V8OSr73lg/4CR9+U2HcZk0UqO8h4GrAUZ9ZP+0rEMQP4K3BPJZ5hjDGpM/VZb3IDzto3iVZ/f7j0fTj0Um/514/Czg2Jf58xOaDCFhxT9VgLjjEJsGEBPH0UlLi6iTod7+w7lSy7NztbOuzZEijrezWc/HDy3mlMFZaoaeLGGJMbSvfBe7/1Jjc16sNp/0rue2s1cjbldPv+RSfZMsbExBIcY4wJ9u0TsGKqt+zEB6FBCuZZ9L0KGrq2oistgfF3J/+9xmQZS3CMMcZt3Rz44j5vWZcToeeFqXl/tRrO4GK3Xz+EZd+Fr2+MCcsSHGOM8dtX4myEua8oUFazodM1Fe+6N/HoMSx0IcHP/mw7kBsTA0twjDHG79cPYNWP3rJTHoF6Kd4rWASG3OstWzEVFo5PbRzGVGEJTXBE5C++T+OKaxtjTIb5+S3v+QGnOa0p6dCuP3Q92Vv2y3tpCcWYqijRLTh3A3cBTRP8XGOMSa49W2FB0IaaR/1farumgvW9yns+d4wzw8sYUyHrojLGGIA5n8A+1+4ojdrFt6FmIrU9ypme7rdrA6yYFrm+MaaMJTjGGAPwyzve8wOHprf1BqBadeg8xFs256P0xGJMFWMJjjHG7NoEC7/wlvUYmp5YgnULGocz52ObTWVMFCzBMcaYOR85C+r5NekMLXqkLx63ToMhryBwvmkRrJ+bvniMqSIswTHGmFlB3VM9MqB7yq9mfWg/wFs29+P0xGJMFWIJjjEmt+3cAIsnessOzJDuKb9up3jP51iCY0xFLMExxuS22e+DuqZeN+8OzbulL55wgtfDWTkdtq1OTyzGVBGW4Bhjctsv73rPM631BqD+ftCqt7ds7ifpicWYKsISHGNM7tq+BpZM8pZlyuypYMGtOJbgGFMuS3CMMblr9vuAa8p1y4OhSce0hVOubqd6zxd9BXu2pScWY6oAS3CMMbkreG+nTG29AWjWFRp3CJyX7oUF49IXjzEZzhIcY0xuKt7p7NDtduBZ6YklGiKh3VQ2m8qYiCzBMcbkpmVTvIv7Ne7o7D+VyYK7qeZ/BiXF4esak+MswTHG5KbgwcXtjkpPHLFo0xdqNw2cF22DpZMi1zcmh1mCY4zJTSEJztHpiSMWefnQ9SRv2a+2+aYx4ViCY4zJPUU7YNUP3rKq0IIDod1Ucz+B0tL0xGJMBkt0glMLqK2q8xL8XGOMSZzlQeNvmnRyFtOrCjoMhII6gfPtq2HVj+mLx5gMldAER1WLVLUokc80xpiEq4rjb/wKakGnQd6yOdZNZUww66IyxuSeqjj+xi24m8qmixsTwhIcY0xuKdoOK6vo+Bu/LkMgr1rgfMNc2DA/ffEYk4EswTHG5JZl33l3D2/SGeq1TF888ajVKDQps1YcYzwswTHG5JYlX3vPq1rrjZ91UxlTrmoVVwklIpcmOpBoqeor6Xq3MSYLZEuC0/Uk+OSPgfMV05zd0ataa5QxSRJXggO8hGcL3pSyBMcYE58922DVDG9ZVRtg7NegNezfyzVFXGHuGOhzWVrDMiZTVLaLSlL0wfVfY4yJz/Kg8TdNu0C9FumLp7K6neI9t+nixpSpTIKTyoTDkhtjTOUtnug9r6rdU37B43AWfeW0Uhlj4u6iOjahURhjTCpU5QX+wmnWDRp3gE2LnPPSvc4O4wednd64jMkAcSU4qvpVogMxxpik2rMNVs/wllXV8Td+Ik4rzuR/B8pmv2cJjjHYNHFjTK5YOhnUtSll065Qt3n64kmU7md6z+d/7mwmakyOswTHGJMbFnzuPW9fxVtv/FodCg0KA+cle2Dep+mLx5gMYQmOMSb7qcK8z7xlnYekJ5ZEE4EDz/CWzX4vLaEYk0kswTHGZL/1c2HrssB5fo2qP/7GrftZ3nPrpjImfQmOiLQSkZ4icqzv01NEWqUrHmNMFpsf1HrT/mioXjs9sSSDdVMZEyJlCY6I5InIUBF5R0TWAMuA6cA432c6sExE1ojIu7661sJkjKm84ASn8wnpiSNZrJvKmBApSSBEZAiwEHgTOANoTuRVi5sDp/vqLhSRLPuXyBiTUnu2wrJvvWWdB6cnlmSybipjPJKe4IjIPcAYoC2BFYnL28fKf01893wiIvclL0JjTFZbNAFKSwLnTTpD4/ZpCydpwnVTzR+bvniMSbN4VzKOiojcDPzZd+pOXLYAPwBzga2+svpAF6A30DCo/p9EZKuqPpTMeI0xWSh49lSXLG0UFoHup8O3TwTKfnkXegxLX0zGpFHSEhwR6QbchzdRmQHcDXyiqiUR7ssHTgHuxEl21HfvvSLyoar+mqyYjTFZprQ0dP2bbOye8jtwqDfB8XdT1aibvpiMSZNkdlHdi5NA+bulHgH6qOoHkZIbAFXdp6ofAH2Bh333K5AP/C2J8Rpjss2ambBjbeC8el0oPDJ98SSbdVMZUyYpCY6I1AJOxklMFHhVVW9Wda+TXj513Aq8TGAA8sm+Zycy1nwR6SEiI0TkcRH5VkR2iYj6Pncn+H0aw+elRL7bmJwTPHuqwzFQrXpaQkkJfzeV2y/vpSUUY9ItWS04RwM1cZKSfcDNlXjWLYC/xaeG79mJ9AbwMzASuBY4HEhoEmWMSZOQ6eFZsnpxeQ4Mmk218EsoKU5PLMakUbLG4LTx/VeB71R1fbwPUtX1IvId0N9XVFhe/TjkB51vAjYCnRP8nmC/EBiAHcmyCq4bYyLZuRFWfO8ty+bxN36tekPdFoGuueLtsGIqtDsqvXEZk2LJSnCauo6XJ+B57mc0ScDz3KYCv+IsNDhdVReLyAicFp1k2qCq7yX5HcbkrgXj8KxI0fIgqL9/2sJJGRHoeBz89L9A2YJxluCYnJOsLqrdruN6CXieewrA7oi14qCq96vqbar6lqouTuSzjTFplIvdU36djveeLxiXnjiMSaNkJTj+aQsC9ErA89zPWBuxljHGgDM9fNGX3rJcSnA6HEtgAiuw5mfYviZt4RiTDslKcKb7/qvAfiJyVnmVyyMiZwLuTTh/qERcxphcsO4X2LUxcF6jPrTqk754Uq1OE2fKuNvCL9ITizFpkpQER1UXAAt8pwI8KSIdY32OiLQHniDQkb5QVecnJsq06yoik0Rko4gUi8haEZkoIneJSMt0B2dMlbboK+95u6MgP6kLt2ce66YyOS6ZC/39k8AifS2Bb0Tk7GhvFpEzgEnAfq7nPJqEONOlJc7MsMZAAc4mo0fjrPS8RESuT8RLRCTsx5istjgowWk/MD1xpFNwgrPwCyjdl55YjEmDZP5K8yxwDXAwTnLSHBgtIrOB14HvgPnANt/1BkAnoB9wHtCDQGKjOGvVPJPEeFNpIfAZMBNnSnot4EDgbKADzno/j4lIY1W9K21RGlMVlRTDkm+8ZR1yMMHZ/1Co2RD2bHHOd2+GVT9C6xzqqjM5LWkJjqqWisgpwGSctWv8e0odCNxTwe3uXccFZz2YU2JZCTmDDVTVieEuiMjtOC04/vVx/iIin6vqpFQFZ0yVt3I67N0ZOK/bApp1S1886ZJfDToe62y46bdgnCU4Jmcks4sKVV2F0w3zJd7WGAhsvxD8IajOV8DRqroymbGmSqTkxndtn6reCTznKr4j+VEZk0XCdU/laresjcMxOSypCQ6Aqq5U1UHApcD3eBOZcPzXvweGA8epaiIWC6xK7iaQ5B1bmf23VDXsx5isFTzAOBe7p/w6DvKer5wOuzalJxZjUixl0wpUdRQwSkQ6A0cCvYFmQENflS3Aepwp5pOzaLZUzFR1lYjMA7rijMdph7PasjGmPEU7nG0J3HJxgLFf/f2gRQ9YO8s5V9/6QD2GpTcuY1Ig5fMmfYnLfJxdwk1krkU8aJS2KIypSpZ9C6UlgfPGHaFhm8j1c0GnQYEEB2DBeEtwTE5IeheViZt7z60t6QrCmCpl0QTveS53T/mFG4dj3dQmB1iCk4FEZD+gi++0GFiaxnCMqTps/ZtQbQ6HgjqB8x1rna0bjMlyluBkprsJDMT+SlV3llPXGAOwc0PQD26B9gPSFk7GqFY99PsQPE7JmCxkCU6CiMhLIqK+z90R6jwgIoXlPCNfRO4BrnYV35fgUI3JTouDVmBoeRDUbpyeWDJN8No3q2akJQxjUimlg4xFpCtwCM7sqfo4WxTERFUrWiQw1pjaA1cEFR/sOj5ORIK/T2+r6o9xvO63wK0iMhn4BpgHbMVZybg7cC7OSsZ+96nqVyFPMcaECu6e6nBMWsLISPv38p6vnpGWMIxJpaQnOCLSCLgVuBhnX6nKSmiCA7Sl/MX0jvZ93BYA8SQ44HQ99fd9ItkF3KqqT8T5DmNyj61/E1lwgrPuV9i7BwpqpiceY1IgqQmOiBwP/BdnRpB7cb94hvD7V0Kuyk7ASWyOwBlE3Azne7MPZ1r4TGA88JKq2mpcxkRr81LYvDhwnlcAhUekL55MU7sxNCiErcuc89ISWPsLtO6d3riMSaKkJTgi0hf4EGehOgjsKwXlr2ScUqo6gQTEo6ojgBEV1PkOZ5NRY0wirZjmPW/TF6rXCV83V+1/SCDBAVj9oyU4JqslswXnGZzkxp/Y7AReA77A6eLZBpREvNsYY6IVPO25lf3gDrFfT/j1w8C5DTQ2WS4pCY6IHIwzmNjfpTQdOE1V1yTjfcaYHOdeqRecGVTGywYamxyTrGnih/n+KzitNGdbcmOMSZo1QQlOix7piSOTRRpobEyWSlaC09z3XwW+VVVbidcYkxw7N8AO1+9P+dWhaef0xZOp/AON/fwDjY3JUslKcHa7jlck6R3GGBM6/qZZN8iPeYmt3LD/Id7z1fGudmFM5ktWgjPPdVw/Se8wxpjQBKflweHrGWegsZsNNDZZLFkJzkTA37lr0xmMMckTMsDYxt9EZAONTQ5JSoKjqjuAl3AGGbcUkdOT8R5jjLEBxjGwgcYmhyRzs807gOU4Sc4TItIqie8yxuSikiLYMNdbZi04kdlAY5NDkpbgqOpmYAiwCmgNfCcipybrfcaYHLR+rvND2q9+a6jVKH3xVAU20NjkiKTuRaWqc0WkN/ACcDLwvogsBMYBc3F20i6N8ZmvJDxQY0zVFDLA2FpvKmQrGpsckfTdxHE2kfwQOBJoCHQCOlbieZbgGGMctoJx7GygsckRyd5NvDvwNs7O2RC6G3i0m1z697Oq6ruJG2MSKbgFxwYYVyzSQOOCmumJx5gkSdoYHBFpD0zASW7cu4i7P1E/LqHBGWOqPlVrwYlHuIHG62ygsck+yWzBeRpoSqDVpRT4HGc38YXYbuLGmMrYtgp2bw6cF9SBRu3TF09Vsv8hsHVZ4HzVj7YDu8k6ydpNvAswmEDX0kLgTFW1XxOMMYkR3HrTojvkJXPliyxiA41NDkjWvwZH+f4rOC03Qy25McYk1JqZ3nPrnope8DicRROgNKYJrcZkvGQlOC18/1Vguqr+XF5lY4yJma1gHL82faFarcD51uWwZGL64jEmCZKV4GxxHS9M0juMMbnMBhjHr0Y96B60g86Po9ITizFJkqwEZ7nruFbEWsYYE4/inbDR/buTQPPuaQunSup1sfd89gfeQdvGVHHJSnAmAUW+Y/u1yhiTWOt+xbMsVuP2UKNu2sKpktoeBQ1d08X3FcGst9MXjzEJlqzdxLcA7+AMMm4vIkcn4z3GmBwVPMDYxt/ELi8Pega14lg3lckiyZxTeTuBsTj/EZGGSXyXMSaXBA8wbnlweuKo6npegGcd1VU/2u7iJmskczfxpcDpOAv6HQBMFJGeyXqfMSaHrPrBe26bbManYSF0OMZb9uNraQnFmEQT1eRs7yQiA3yHXYCHcDbaVJztGz4D5hDfbuI2l7ECIlL2h5qsP19j0mbnBni4E54xODfNg3otIt5iyvHzW/D2FYHz2k3gxjlQrXr6YjImSiKBFkhV9WzrlMytGiYQfnPNY3yfeCip2QHdGJOpFn6J55+WlgdZclMZ3U6Bmg1gz1bnfNdGmPdp6DRyY6qYVKxr7s+oFG/CE7zxZrQfY0wuWzjee97p+PTEkS0KasFB53jLZlg3lan6kp3gRNpF3BIVY0zsSkthgSU4CRe8Js78z2DryvTEYkyCJLO757IkPtsYk4vW/gw71wXOq9eD1n3TF0+22K8nND8Q1vlmUGkpfPMvOPmhtIZlTGUkLcFR1ZeT9WxjTI4Kbr1pP8AGwyaCCPS9Cj66IVA2/SU4+kao1zJdURlTKakYg2OMidMtt9yCiCAiXHzxxRXfALRr167sHvcnPz+fRo0a0bdvX+666y42bNgQ8RkjRowI+4zgT7t27cqN5bHHHkNEuOSSS8rKjjnmmLL7q1WrxsqV5XeFFBUV0aRJE+eeATfR7rHtgYudBgGwZMmSiDFWr16dli1bMnjwYJ555hn27t1b8TfR56mnnip7zlFHHRXVPe6vz/2pU6cOHTt25Pzzz2fs2LFRx5AyPS+E+q0D5/uKnFYcY6ooS3CMyVAlJSW8+uqrZedvv/02W7Zsifr+mjVr0qJFi7JPvXr12LJlC9OmTeOee+6hW7duTJ8+vdxn5OXleZ4R/GnWrFm597///vsAnHHGGWGv79u3j1deeaXcZ7z33nts2rQp/EVfguNWv359T4zVqlVj7dq1jBs3jt/85jccfvjhkZ8X5MUXXyw7/uabb5gzZ05U9wEUFBR44iguLmbRokWMHj2aE088kauuuiqzlnGoVgOOusFb9v2LsH1tWsIxprIswTEmQ3388cesWbOG7t27c9xxx7Fnzx7++9//Rn3/eeedx5o1a8o+W7ZsYdu2bTzzzDPUrFmTjRs3ctFFF5X7Q7ZNmzaeZwR/pk2bFvHeTZs2MWnSJGrUqMEJJ5wQct3f+jNy5Mhyvw7/9XatgqaCN+kMjdqF1P/Xv/7liXHXrl2sWLGCG264AYAffviB2267rdx3Avz0009Mnz6dRo0acdFFFwHehKciRx55pCeOPXv28P3333P00c7ONc8//3xMz0uJQy+FevsHzkv2wOR/py8eYyrBEhxjMtQLL7wAwKWXXsqll14KxPYDNpx69epx9dVXc8cddwAwd+7cmFolYvHxxx9TUlLCcccdR7169UKuDxgwgHbt2jF//nwmTZoU9hkrVqzg888/p27dugw7rJX3YpjWm0hatWrFo48+yqBBzj3vvfdehff4v//nnXceV111FQCvvPIKJSUlUb/XLT8/n969e/P+++/TpEkTzzsyRrhWnGkvwI71aQnHmMqwBMeYJFmzZg1/+tOfOOSQQ2jQoAE1a9akQ4cOXHnllcyePbvce1evXs2YMWPIy8vj4osvZtiwYdSpU4fp06fz008/VTq2nj17lh3v2LGj0s8Lx59EROqeEhFGjBgBRG7FeemllygtLeWcs8+m7u7l3otxTA/3f90Vfc1FRUW89pqzFszw4cMZMGAA7du3Z+3atXz88ccxv9etUaNG9OvXD4BffsnAfZ8OHQ51XQOLS3bDt4+nLx5j4mQJjjFJ8NFHH9G5c2f+/ve/M3PmTHbv3k21atVYvHgxL7zwAr169Sp37MnLL79MSUkJgwYNolWrVk4LxrBhQGJ+6/cnSfn5+XTs2LHSzwu2Z88exo4di4hw+umRV8T1D2Z+44032LlzZ8j1l156CYDLhg4JrLQLkF8D2vaPOS7/1921a9dy67377rts2rSJLl26cPjhh3sGSifi++/vFty3b1+ln5VwBTWh//XesqnPw86N6YnHmDjZtgcmau3+VLnfXDPVkgdPSejzpk6dyrBhwyguLuaaa67hhhtuoHPnzuTn57Ns2TL+/ve/89RTT3HFFVfQvXt3+vTpE/IMf1eUv2sKnJaEV155hddee42HH36YGjVqxBzbjh07GD16NPfddx8AV1xxBY0bN47zK41s/Pjx7Ny5k379+rHffvtFrNe2bVuOO+44xo8fz5tvvlnWogPw1VdfsXDhQjp16sTRzbbimSDe9kioXjvqeFavXs0jjzzCuHHjALj55pvLre9PYtyzvy699FLuuecexowZw+rVq8v9usqzefNmpk6dCkCHDh3iekbS9R4Bkx4NrDm0dydMeQoG3ZnWsIyJhbXgGJNg1157LcXFxdx55508/fTTdOvWjfz8fAAKCwt58sknue666ygpKeHee+8NuX/ixInMnz+funXrMnTo0LLyY489ljZt2rBp06aoxpCMHj2ali1bln0aNWpEvXr1uPLKKyksLOT+++/nqaeeKvcZy5cv9zwj+LNt27aw91U0e8rtssucNUGDu6n8Sd5ll10GC8Z5byqne+r666/3xFinTh32339/HnnkEY4++mjefvttLrjggoj3L126lPHjx4dMb+/YsSP9+/enpKSkwplf4ezbt4/p06dzxhlnsHHjxsDXlomq14b+13nLfv0wPbEYEydLcIxJoJ9++olp06ZRUFDATTfdFLGev2Vm3LhxId0U/taDs88+m9q1A60UsXaT7Nmzh7Vr15Z93FPMN23axPr169m1a1e5zygtLfU8I/hTWloaco+q8uGHzg/DaBKcoUOH0qBBAyZOnMjChQsB2L59O2+//TZ5eXkMv/A8WBI0CLmcBGfbtm2eGN1f4/r161m5cmW5M8defPFFVJWBAwfStm1bz7Xhw4eX1anI5MmTPYlWzZo16dOnD19//TUAw4YN49prr63wOWlz6KUgrh8RG+bC9jXpi8eYGFmCY0wC+WcDlZaW0rVr14gtHyeeeCIAO3fuLPttHpwfzm+99Rbg7Z7y8/+AHT9+PEuXLi03luHDh6OqZZ/S0lI2btzIBx98QPPmzXn00Ufp378/69dHniHTtm1bzzOCPw0bNgy5Z8qUKaxZs4ZOnTrRvXv38r9hQK1atTj//POBQCvO6NGj2blzJ0OGDKGVrnIGuvrlVYNmkcfQjBw50hNjSUkJK1as4Omnn2bt2rVcd911XHjhhWGTnNLS0rJxP+G+/+eeey41a9Zk3rx5EWd++e3du9eTaPlnX4kITz75JG+99RYFBQXlPiOtajaA/Xt5yxZ/nZ5YjIlHef942adqfgjs3K4mtW6//XZ1f/+j+SxevLjs/meeeUYBLSws1NLS0rDvOPzwwxXQu+++O+z1tm3bKqDDhw+PGOfy5cu1fv36Cuhll10Wcn348OEKaNu2bWP58lVV9dZbb1VA//jHP4a9PnDgwJD4vvvuOwW0TZs2um/fPj3yyCMV0DfeeEN14iOqd9XXuwZWd2JqVjfkmYsXLy77fo4cOTJibB988EFZvZdffjnk+tixYxXQ2rVr67Zt28I+4/zzz1dAR4wYUe7XN3DgwLKy4uJinT9/vt58880qIlqnTh2dMGFCxDgzxud3q95VP/B57/fpjsgYj6Cfd56fhdaCY0wC+bubunXrFnVC6t7uwN/1tGzZMvLy8sIu+T9lyhQg0FIRj9atW5e1Ir3xxhtxr+0STizjb/z69u1L9+7dWb58OU8++SSTJ0+mcePGzgysFd97K1eLfXC132mnnUbz5s0Bwi6a6P/+79q1i/r164f9/r/++usAvPnmm2zfvj3kGeEUFBTQqVMnHnroIe666y527tzJueeey7p16yq+OZ3aD/CeL56YnjiMiYMlOMYkUMuWzvohixYtCjvtuTyzZs0qm10TDf9g2Hj5x5fs3Lmz3G6qWMybN485c+bQrFkzjjzyyJju9Q+4/eMf/wjAhRdeSI3q1WFF0GrJ+fEnOBD4uhcvXuwp37hxY1lyFo2dO3cyevTomN9/++2307FjR9atW8edd2b4rKQ2/SDftZnplqWweUnawjEmFpbgGJNA/fs7a7MUFxfz7rvvxnSvv/Xg0EMPZfv27eV+zjzzTM898VixYkXZcZ06deJ+jpt/dtepp55KXl5s/7xccsklVKtWjeLiYgAuv/xy2Lo8MFXZL79y41b8X3fw1zxq1CiKiopo3rw5W7duLff7f/31zjox8Xz/CwoK+POf/1x2/7x58yr19SRV9drQuq+3zMbhmCrCEhxjEqhPnz706uUMzLzjjjsqbBnxb/pYXFzMqFGjAGcga926dcv9nHfeeYCzIN3mzZtjjnPDhg18/vnnAHTu3Jn69evH/Ixw/C0g/gQsFi1atODRRx/lpptu4q9//avzfQxuvamkCRMmsHr1aoCQ9Yf8M6OGDh1K/fr1y/3++wdFT5kypcJVqcO5+OKLadu2Lfv27eOvf/1rJb+qJAvppvoqPXEYEyNLcIxJIBHh6aefpkaNGixbtox+/frx1ltveaYqr1y5klGjRjF48GBuvfVWwEkMNmzYAMA555xT4XtOO+00atWq5dlSIBqlpaVMmzaNM888s+x9N954YyxfYkTr1q1jypQp1KpVi8GDB8f1jGuvvZZ//OMf/OUvf3EKgsffxKmoqIj333+/bNPMatWqcd11gXVepk2bxsyZMwEnwazI4YcfTmFhIRDf/mDVqlUrW2zw9ddfjytJSplw43DiHPtlTCpZgmNMgvXt25cPP/yQJk2asHjxYs455xzq169P06ZNqVOnDq1bt+aSSy4pW1UXAl0dvXv3jmp12zp16nDSSSd57g0WvNBfixYtqFOnDn379uWbb74BnOTmN7/5TWW/ZAA+/PBDSktLGTJkCLVq1UrIM+NpwQle6K958+bUrVuXM888k1WrVlG7dm1GjRpFjx49yu7xfw+bN2/OgAEDIj3a4+yzzwacDTj37t0bc5xXXHEFLVu2pLS0lLvuuivm+1OmVW8ocHXn7VgLGzK4W80YH0twjEmCwYMHs2DBAh544AGOOuooGjRowJYtW8jLy6N79+5cccUVfPDBBzz++OMsX768rLsomtYDP3/dGTNm8MMPP4RcD17ob926deTl5dGlSxcuv/xyJk+ezCOPPJKYL5j4Zk+Vq6QIVs+M+bbghf42bNhA7dq16d27N7fccgu//PJLWRcfwO7du8tmRg0bNqxs1emK+L//69evL1vYMBY1a9Ysaz17++23E7KJalJUqw5tj/CW2WwqUwVIvNNMRSTVu8QpsAPYDqwEfgG+Bd5W1dgHIWQxESn7Q433z9eYWOzatYumTZtSVFTEmjVraNasWeUfumI6PH9c4Lx+a7gxA3ffzgXf/As+/0vg/IDT4LxR6YvHGB8RKTtWVXFfq0wLjqT4kwfUB1oBfYERwDPAGhF5UkQaxPVFiOSLSA8RGSEij4vItyKyS0TU97k7nudG8d46InKjiHwjIutEZI+ILBWR0SJyQjLeaUyyjB07lt27d9O/f//EJDcQ2j3VOnRTUpMiIeNwvoYw23QYk0kqu5t4upsHBCgAfgMMEpEBqhrryllvAEMrrJVAItILeAsIHmxR6PucKyL/BS5T1eJUxmZMPOrUqcNdd93F0UcfnbiHhiQ4hyXu2SY2LQ92tm7Ys9U537MF1v4M+x2S1rCMKU9lExypuErSKU4cXYDRwLEx3h/c4b4J2Ah0rnxooUSkLTAGaOErmgqMAjYABwFXA02AC4FS4JIwjzEmowwZMoQhQ4Yk9qErg2ZQWQtO+uTlQ7ujYc5HgbLFEy3BMRkt7i4qVc1L9QeohfPDvydwKeBfSc3fkjRARGL9V3Yq8CBwDtBBVZsA98f7fYnCYwSSmxeBI1T1cVX9n6reDhwKLPNdv1hETkliLMZkph3rvSvm5lWzH6bp1n6g99wGGpsMV9kWnJRS1SKgCNgMzARGichg4EOcriqAYcBnMTwzmcmMh4gcApzpO10G/F5VPR3ZqrpMRH4LfOwrutt1bExuCG69aXkQFCRo6rmJT/A4nKWTYd/eSq8sbUyyVPlp4qr6OfA0ge6y/mkMpyLnuY6fVdU9EeqNARb4jvuISMfkhmVMhgkef9PKuqfSrllXqNM8cF68A1aGLk9gTKao8gmOz3u+/wqB7p9M5O4++zRSJd8W8GNdRTaryuSW4BWMbYBx+onYtg2mSklagiMi7ZPwzHsjXPLvGqhAw0S/NxFEJA84wHdaAlS0qpf7X/geEWtV/N6wH2MyVum+0JYBG2CcGToEjcNZZAmOyVzJbMH5VEQaJ+phIvIX4LYIl/3rpPvXy8lErYHavuOVqlpSQf2lruMuyQnJmAy0fi4Ubw+c12oMjSvevsKkQPBA4xVToXhX+LrGpFkyk4HOwMciUrOyDxKRm3EG24alqktds62iW2c99Rq6jjdEUX9jhHuNyW7hpodbq2NmaNQWGrYNnO8rhmXfpi8eY8qR7NaOvsBoqUSfiIhchzONu6qr6zqONLjYbbfruF6CYzEmc9kCf5ktuJvKpoubDJWK7pxTgf/Ec6OIXIOzbozf1kQElAFStgK0qob9GJOxVs3wnrfqnZYwTAQh6+HYOByTmZKZ4PyHwD5SV4nIn2O5WUSGA08RWKl4O3BSooNMoR2u42gW9HDX2R6xljHZZN9eWD/HW2YL/GWW4ARn1QzYbfsdm8yTzATnWuAD37EAfxWREdHcKCLnAy+47t0JnKqq3yU6yBTa4jpuEkV9d50tkSoZk1U2zHfGdfjVbQl1mqYvHhOqbjNofqCrQGHJpLSFY0wkSUtwfGu5nA/4kxIBnqlop2wRGQq8QqD1Zzdwpqp+naxYU2QF4J9u0FpEKlpF2jWSj3nJCcmYDLN2lve8ZdwrJJhksunipgpI6hgc30q9p+Ksyqs42ym86dtNO4SInAr8D2cDTMHZluFsVR2fzDhTwbclw6++02pARe3u7oU/ZkWsZUw2WfOz97yFJTgZyRb8M1VA0gcZq+pG4ERgPU6SUxf4RETauev5Nsl8EycJEpy1bS5Q1THJjjGFolqd2DfrzH19bKS6xmSVkBacg9IThylf2/4grhU5NsyDbavTF48xYaRkUTxVXYTTkrMLJ8lpAYzxLwQoIsfg7Axew3e9BLhUVd9LRXwp9Ibr+Jpy1gg6CejkO/5eVRcmNyyTqW655Zay1acvvvjiqO5p165d2NWr8/PzadSoEX379uWuu+5iw4bIyzGNGDEi4irY7k+7du3KjeWxxx5DRLjkkkvCXt+2bRuPP/44p512Gm3btqXOFe9T675ttP7ndk5+bRf/fPd71qxZ47nn7rvvjiq2WGMF6Nu3b1n9559/vsL6S5Ysifi9btiwIX369OHWW29l2bJlFT6rSqlZH1od6i2z6eIm00SaRpyMD3AyUAzs832+AQbhzBIq9X1KgEtSGVeYOEfgJFoK3B3lPS9Fcw9OIuev9zyQF3S9EGcVY3+dU+KI33+vmqpr79692rJly7I/y5o1a+rmzZsrvK9t27Zl9Vu0aFH2adCggbr/32jSpIl+//33YZ8xfPhwBTQvL8/zjOBPnz59yo3lmGOOUUDffPPNkGvPP/+8NmrUyBNTrWpogxp4ymrWrKn33ntv2X0PP/xwuTG5P+7nd+3atdxYf/75Z897Dz/88Aq/14sXLy6rX79+/bL3Nm7c2POsOnXq6CeffFLh86qUcfeo3lU/8Hn3t+mOyOSgoJ933p+FwQXJ/gBX+hKZfa5PqavsqhTH0x64N+jzgeubNjHM9V5hnhNtgtMWWOOqOwX4Pc6A7PtwVjn2XxsV59dkCU4WeO+99xTQ7t2763HHHaeAPvnkkxXe509whg8fHnJt27Zt+swzz2jNmjXLfuiXlpaG1PMnOG3bto07/o0bN2q1atW0Ro0aum3bNs+122+/vez/0QMOOEBffvllXTvlzbIflrtur6djr+uhF198sebn5+shhxwS8/tLS0v1jDPOUEDz8/P1s88+K7f+DTfcoICOGDFC69Wrp4D+8ssv5d7jTnBGjhzpubZz504dOXKkNmzYUAFt0KCBbty4MeavI2Mt+sqb4DzSXTXM/0vGJFNGJThOPPzVldCUuj5/SEMsx7i/QVF+RoR5TlQJjq9uL2BhBe94Dage59dkCU4WOO200xTQBx98UF966SUFtHfv3hXeV16C4/e3v/2t7P+R2bNnh1xPRILzyiuvKKAnnXSSp/z1118ve/c555yje/bscS5Mesz7A/O936mq07Jy4YUXxvz+u+66q+w9//jHP8qtW1RUpE2bNlVAp0yZoiNGjFBAb7rppnLvKy/B8fN/HwB95plnYv46MlbxbtW/Nff+mW1YkO6oTI4pL8FJy8aUqnoXTkLg38JBgVtU9fF0xJNqqvojcDBwE/AtTqtNEbAcZ6D1Sap6kaoWR36KyXRr1qzhT3/6E4cccggNGjSgZs2adOjQgSuvvJLZs2eXe+/q1asZM2YMeXl5XHzxxQwbNow6deowffp0fvqpoo3oK9azZ8+y4x07dkSuWAnvvfceAGeccUZZWXFxMTfffDMA3bt355VXXqFGjRrOxTVBA4xbOAOMe/TowauvvhrTu99//33uueceAC688EJuuummCutv2LCBrl270q9fP4YPHw7AK6+8wt69e8u9tyInnnhi2fEvv/xSqWdllIKa0Kaft2zRhLSEYkw4Fa3FEpaIDKi4VoVeAwYC7YBPgKnRPFdVEzqSTVUnEEi0KvOcEThjd6KtvxP4p+9jssxHH33EBRdcUJY8FBQUUL16dRYvXswLL7zAq6++ynPPPcell14a9v6XX36ZkpISBg8eTKtWrQAYNmwYr7zyCi+88AL//ve/KxWfP0nKz8+nY8eOlXpWOHv27GHs2LGICKeffnpZ+fvvv8/y5csBuO2226hZ0zXOvpw1cPLyov9dbM6cOVxyySWoKr169YpqsPALLzjrivr/PAYOHEi7du1YsmQJH374IUOHDo36/eXZt29fQp6TMToM9E4RX/wVHHZF+uIxxiWuBAeYAAnbT0lwBh+fHEVdJf6YTWXd3SDdESTH3Ynd4mzq1KkMGzaM4uJirrnmGm644QY6d+5Mfn4+y5Yt4+9//ztPPfUUV1xxBd27d6dPnz4hz3jxxRcBPAnQ8OHDeeWVV3jttdd4+OGHAy0fMdixYwejR4/mvvvuA+CKK66gcePGcX6lkY0fP56dO3fSr18/9ttvP085OAmLO/GhpMiZauzW4kBitXXrVs444wy2b99Os2bNeO+996hVq/ydUZYvX87nn3/umanmn/n1t7/9jRdffLFSCc6YMYGVLjp06BD3czJS+2OAewLni7+G0lKIISE1Jlkq+3+hVPIDgUQplnuMyVjXXnstxcXF3HnnnTz99NN069aN/HxnzZDCwkKefPJJrrvuOkpKSrj33ntD7p84cSLz58+nbt26nh+sxx57LG3atGHTpk1l3T/lGT16NC1btiz7NGrUiHr16nHllVdSWFjI/fffz1NPPVXuM5YvX+55RvBn27ZtYe97//33AW/3FAS6aDp27Ej9+vUDF9bPgdKSwHmDNlCrUYVfo5uqctFFFzFv3jyqVavGm2++SWFhYYX3jRw5ktLSUo499lhPfX9y+emnn7Jy5cqYYgHYtWsXL7/8Mtdffz0ANWrU4IILLoj5ORlt/55Qw/WLz+5NsPbniNWNSaXKJDiJSDYscTFZ5aeffmLatGkUFBSUO+7D/8Nz3LhxId0W/u6Ss88+m9q1a5eVu9eT8dcpz549e1i7dm3ZZ8uWLWXXNm3axPr169m1a1fkBwClpaWeZwR/SktLQ+5RVT788EMgNMHZuHEjQGirUcj4m9hXML7zzjv5+OOPAfjnP//JwIEDK7xHVRk5ciRASHdhp06dOPLII9m3bx8vv/xyhc+6/vrryxK/Jk2aUKdOHUaMGMGWLVsoKCjg5Zdf9rRmZYW8fGh3lLfMtm0wGSLe7p7LEhqFMVli0iRn08HS0lK6du0asZ4/qdm5cycbN26kefPmgLPw3VtvvQWE/sAFp5vq/vvvZ/z48SxdupS2bdtGfMfw4cN56aWXys5Vlc2bN/PNN99w22238eijjzJu3DjGjx9Ps2bNwj6jbdu2LFmypNyvOdiUKVNYs2YNnTp1onv37p5rzqQHJ1nzqOQeVO+88w73338/4CxS+Ic//CGq+8aPH8+SJUuoU6cOw4YNC7k+fPhwJk+ezMiRI7nttttC43bZtm1b2BatwsJCxo4dS7du3aL8aqqYDgNh7seB88VfQf/r0hePMT5xJTiqWvGvMyb7JHisSjZatWoV4CQwa9eujeoedyvK66+/zq5duygsLOSYY44JqdulSxcOP/xwpkyZwksvvcRdd90VdWwiQuPGjTnttNPo1asXBx54ID///DO33npr2ZifRPB3T5155pkh15o2dXYG97fklKnEHlS//PILw4cPR1Xp27cvTz/9dNT3+r/us846i7p164ZcP++887j++utZsGABEydOLLdVaOTIkYwYMQJwkp0ffviBP//5z3zzzTdcdtllfP7552HfUeW1D/qeLJ0MJcVQrXp64jHGx0aCGZNA/paZbt26Rb1ukXsLAX/X07Jly8jLywu7DcCUKVMA5weqv0UkVq1bty6bvvzGG29QUlJSwR3RizT+BuDAA52BwwsXLgy0dqjGvQfVli1bOPPMM9mxYwctW7bknXfeiXrw9ebNm3n33XcBGDVqVNjvdcOGDdmzZw8QXbegX/369TnmmGP47LPPOPDAA5kyZQrXXntt1PdXKc26Qt2WgfO9u2Dl9+mLxxgfS3CMSaCWLZ1/6BctWsTOnTtjunfWrFlMnTo16vpLly4tm5UUD3/31s6dO1m/fn3cz3GbN28ec+bMoVmzZhx55JEh1wcNGgQ4XXgffPCBU7htFezeHKhUUAcata/wXaWlpVxwwQUsWLCA6tWr8/bbb5dNqY/Ga6+9Vpa8ROPtt9+OOKg6ktq1a/P4487yXi+//DKTJ0+O6f4qQSR0d3Ebh2MygCU4xiRQ//79AWdBO3/rQLT8LQSHHnoo27dvL/fj7/6JpVUh2IoVK8qO69SpE/dz3Pyzu0499dSwa9ecccYZtG7dGoAHHnjASTCCW29adA+ZZhxuMPPtt9/Op59+CsDjjz8eNqEqj/97d/3115f7vd66dSvNmjVj165d/O9//4vpHeDMfvN3bd16660x318ldAjqplpsCY5JP0twjEmgPn360KtXLwDuuOOOCltGNm3aBDgJ0ahRowA499xzqVu3brmf8847D4B3332XzZs3R3x+JBs2bODzzz8HoHPnzt4p25VQ3vgbgOrVq/PQQw8BMHv2bIYPH07x8h+9lYLG38yePTtkwPWbb77J3//+dwCuueYarr766pji/OGHH5gxYwYAF1xwQbnf6/r165dN1483obzjjjsAZxC6//ueVYLH4ayYBkXJWSHbmGhZgmNMAokITz/9NDVq1GDZsmX069ePt956yzOQeOXKlYwaNYrBgweX/Ubv3yoA4JxzzqnwPaeddhq1atWiqKiI1157Ler4SktLmTZtGmeeeWbZ+2688cZYvsSI1q1bx5QpU6hVqxaDBw+OWO+CCy7glltuAZzxP71GPMirPxWzfqevlaZlD/bs2cP48eO5/PLLOeSQQ5g1K9DKM3PmTC67zJnI2b9//7IuoFj4E5XCwkL69etXQW0n6QSYNm2aJ5ZoDR48mMMOOwxwprNnnYZtoLFrRezSElj2bfriMQbSs9mmfZL7wTbbTLvPPvtMmzRpUvbnkJ+fr02aNNHatWt7Noe78sorVVX1hBNOiHozTb+hQ4cqoD179vSU+zfbrFmzprZo0aLs07x587JdxP2fG2+8Meyz49ls8/nnn1dAzzjjjKjqP/PMM2U7bfs/tQvQhvXrectq19aHHnqo7D7/9wrfDt3ur7Giz7Jly3T37t1l761oM02/kpISbd68uQJ6ww03lJVHs9mm37vvvltW96OPPorqvVXKhzd4N9789PZ0R2RyQNDPO+/PwuAC+1T9jyU4mWHz5s36wAMP6FFHHaWNGzfW/Px8rVu3rnbv3l2vuOIK/eCDD3T37t26bNkyzcvLU0D//ve/R/18967c06dPLyv3JzjhPrVr19YuXbro5ZdfrpMnT4747HgSHP/u5y+++GLU92xZt0ofO7GWnty5mrapL1qrGlqjRg1t1aqVnnTSSfqvf/1L169f77ln4MCBEb++ij6LFy/WUaNGlZ1/9913Ucf6m9/8RgFt2rSpFhUVqWpsCU5paan26NFDAT300EOjfm+VMetdb4Lzn/7pjsjkgPISHHGum2wiImV/qPbna1Jh165dNG3alKKiItasWRNx4cAQy6fCC67urEbt4foZSYnRJNmuTfBQBzzbFN68COo0SVtIJvu5F99UVc9KnDYGxxhTaWPHjmX37t30798/+uQGYEXQein790xoXCaFajcOXb9oycT0xGIMtjO3MSYB6tSpw1133cXRRx8d243BC8K1PixxQZnU6zAQ1swMnC+aAAeelbZwTG6zBMcYU2lDhgxhyJAhsd8Y3ILTqk9iAjLp0f4YmOya1WYL/pk0si4qY0x67FgPW5YGzvOqwX4Hpy8eU3ltj4C8gsD55sWwembk+sYkkSU4xpj0CO6eankQFNRKTywmMarXcZIct6/+np5YTM6zBMcYkx7WPZWdDv+993zOR9aKY9Ii7gRHRIaKSGI2sDHG5J4V07znrS3ByQpdToD9e3nLJjyYnlhMTqtMC85bwAYRGSMivxOR1okKyhiT5UpLYVXQHlQ2gyo7iMAxt3nL5n4Mq2akJRyTuyrbRVUdGAI8DiwVkR9E5K8i0rvyoRljstaGeVC0LXBeqxE07pC+eExidR4CrYJ+DFgrjkmxyiQ4uwEJ+hwC/BmYKiIrRORpETlZRGpUPlRjTNYI7p5q1dv5zd9kh3CtOPPGwMof0hOPyUmVSXCaAKcDzwGrfWXuZGd/4CrgQ5yurHdF5HIRaV6JdxpjsoEt8Jf9Oh0fOnDcWnFMCsWd4KjqHlX9SFWvUdVWQD/gXsA9XN6f7NQhkAytEpHJInKbiBxYidiNMVXViunec5tBlX1E4NigVpz5Y0P/7I1JkoRNE1fVaar6F1XtCbQF/gB8DhT7qviTnTxcyZCILBSRx0TkOBHJT1Q8xpgMVbQD1v3iLWt1aHpiMcnVcRC07ust+/7F9MRick5S1sFR1eWq+qSqngA0Bc4FXgU2+qq4u7LaE0iGNojI/0TkAhFpmIzYjDFptnoGaGngvEknZ6NGk31EYMAfvWW2AWdqlJbCl/fDowfB6xfB7i3pjijlkr7Qn6ruUNW3VHU40AIYCPwDmOeq5k92GuAkQ6OAdSLyhYj8n4h0THacxpgUCRlgbN1TWa39AMivHjjfsgw2L41c31SeKnx6q7OK9NZlzmKLk/+d7qhSLqUrGatqqap+raq3qGo3oCtwM/A1sM9XzZ/sVMOVDInILyLyoIgcKWLTLYypsoJXMLYF/rJbQa3QQeRLv0lPLLliwoMw9Vlv2bzPUhvDrk0w8WH4+p9QvDO17/ZJ61YNqjpfVR9R1YFAc+BSnAUEt/uquLuyuhFIhtaIyIsi0ivMY40xmUrVEpxc1O4o7/mSSemJIxdMeRq+CjNbbe2s1HZTvXUZfHEvjP8rvDoU9pWk7t0+GbMXlapuVtVRqnouzridE4AnAH9bpj/RAWgGDAdOS3mgxpj4bVsJO9YEzqvVhBY90hePSY2QBOfr0Do7N8Cst2HjwtTElI1+et3pmgpLYfnU1MSxfi4smhA4Xz4FJv8rNe92yZgEx01V96rq56p6naq2B3oCdwLTyr/TGJPRgltv9jsE8gvSE4tJndaHhY7D2bIscL5jPTx9FLx1OTx1ROp+EGeTxRPhvd+VX2fZ5MS8a+8e+O5Z+Pap8K1Cv34YWvblA7BmVmLeH6WMTHCCqepMVb1PVfvhLCB4DfAxsCu9kRljYmIL/OWmglqhg8mXuMbhfPc0bPetF7uvCD67M3WxZYsv7wfdFzjPK4BDh3vrLP228u9RhdEXw5ibYextTleUqrdOuASndC+8+xsoKQ69liRVIsFxU9W1qvqcqp6uqv9IdzzGmBgsDfoNMni/IpO9Io3D2bcXfhzlvbZ8Ciybkpq4ssHePaGto8Oeg2P+5C1bOR327q7cu35+ExZ8Hjhf+IX33VuWOUtBhLP2Z5j4UOXeH4Mql+CYKmzbKpjzCexYl+5ITDrs2Rq6g3jb/umJxaRepHE488Z6x2X5fZN705rjtnaW00Li17AQDjwL6u8PDdsGykv3OklOvPZshbF3hJZ//0LgOFzrjdvX/0zZataW4JjkWz8P3v0tPHYQvH4BPN7HKTO5Zelk7wJ/zbpBvRbpi8ekVsg4nKXOb/vTXwpff+7H9u9EtIKTFnfLaNsjvdcq00315f2wM8wvqLPecaaFQ2iCM+AWqNsycK774L3fVL4lKQqW4JjkWTUDRl8CT/aFn/4Lpb5pgkVbYdxdaQ3NpMHioBVs2w9ITxwmParXDu2S/Ol1WDAu8j05uDhdXIITnP1dW58UHuG9Fu9A49U/ha6t47evyOlm3L42tGvx0Evg9Me9ZRvmOVPIk8wSHJMcP74Gzx4Dv34AaOj1uZ84CZDJHZbgmOBuqokP4/n3oXpd7/WZo2Hb6qSHVeXF0oKzfGrsa9KUlsLHN3lbYAlab/f7F2HOh3j+PPfr6XSXdRkCh17qrT9ztNPllUSW4JjE274WPr6RsImN21epG2xm0mznBmecQBmx8Te5KDjB2Rc0o2bwX6FRO+/1755OelhV2u4tsHFB4FzynOUX/Jp0gjrNAufFO5zBvrGYMSp0i5UznvR2OW5e7Kyg7HaAa6m6IfdBg0LnuNup8NtvoWaD2OKIkSU4JvEm/xtK9njLWh4ER93oLZv7sdPsabJf8MJu+x1iG2zmotZ9nenL4RTUgYPPgyP/4C3//kXYsy35sVVVwQP3mx0ANVwtYSJQeLi3TizjcHZtgs+DhhR0PQV6XQTdz/CW71zvPT/g9MBxzfpw1n/grGfgvFFQtxnJZgmOSaydG5x/kNyOuxOu+RoG/SW0D37C31MXm0kf654y4IzDibQ1x0FnQ4160PMiqN00UF60DaaPTE18VVFI99ShoXUKgwcax7AX2NePwO5NgfNqteAkX0tNnysi39e0KzTr4i1rdxQccr6TdKWAJTgmsb59Ava61l+s2xKOuNb5H1oEjrnNW99acXJDSIIzMD1xmPQL7qby6z3C+W9BLeh3jffalP9A6b6QW7LWjnUw47+w8MuKx8us/MF7Hi7BaRs80HhK6OJ84WxdAVOf85YN+KMzrgaclqHm3cPf2/308OUpZAmOSZxdm0L/MvS/DgpqBs47HR/aimNjcbLb1pXeMQJ51UKbzE3uCDf2quXBsL9r7+TDroSC2oHz7atDx4Bkoy3LnMG8j/aA934Lr57pLK8x/m+waXH4e8obYOzX4iDvAO5dG2DD/IrjmfCgM0PKr34rOOL3gXMROCxCK84B6d8q0hIckzhT/uMMYPOr0wx6X+atIwIDg1bXnPORteJks+DWm1Z9vGMETG5pE2YcTu8R3m6L2o2h68neOnPHJD20tNmwwFkr7N+9YNrz3qRi+yr4+h/w757w0qmw3JXobVsVunltuBaV/GrO992touni6+fBjNe8ZQNvdVrY3A4+L3T2W8NCJ2lNM0twTGLs3hI62+HIPzh97sE6D/au0wDWipPNbPyNcateB7qeFDiv3RQOOie0nrsOOCseZ6OFXzobjbrXCotkydfw6lmw3ZfUBLfelLd5bawL/n15r3daeJNOzvioYDXqwcHnessOOD1l42zKYwmOSYypzzqDAf1qNY48AC3cWJx5n0LxzuTFZ9JD1RIcE+rkh6HXxc504QvfcGbYBOs0CCQ/cL7+V9i8JGUhpkTxTmcH8JIwq/rW28+ZWRZyz3aY9JhzXN4Cf8GCBxov+NxpFSstDa27cjrMft9bdtydTktQOANuDqxWXLtp6Ey4NLEEx1Tenm3w7ZPesiN+X343ROfB3vUuSkts4b9stGkRbFsROK9W03YQN1CvpbOOyvmvQesIG67WahS6Cu/cT5MfWyp9/U+nC8qtSSc44ym4fib8cS6c9m/v+CRwZqpuWxXd+Bv3Nfe6Nbs2wv/Oh6f7w8w3vFsnjL/He+9+PUOnhLvV3x+unQbDP4Lrf3L+fDOAJTim8n79EPZsCZzXbAB9ry7/HhFoEzTQdMXUhIdm0iy49aZNP++gc2PK0/VE7/m8LBqHs2lR6FYUfa+G30911pipVt3p/uk9HC4f6wzw9dtX5EzfDv6lMNwMKr+CmtBxUGj5utnwzlVwX0u4bz/4Z3dYNMFb5/i7K+5yqlkf2h+dUePrLMHxEcd5IvKRiKwQkSIRWS0i40XkShGJ0DYX17s0hs9LiXpv0gTPbuh9Wfgm52DB62Gs+D5xMZnMYN1TpjK6BI3DWfJN9iz6N/bP3pWc67Zw1grLyw+tW60GHH2Tt2zaC95hATUbQuMO5b/ztMegy4mRr+/dBdtWesvaD4SOx5b/3AyV8gRHRKqJSHMRaS8iXUSktYgkd73mimNqBIwDXgdOAVoB1YGWwHHAc8B3IlKYtiAzWfBKmpHWuQgWPKp/+dTo1mYwVcPGhbDoS2+ZrX9jYtG0EzTuGDgv3QsLv0hfPImyYJyzBpjb4HucFptIel0CDdq4CoL+rWzVu+JWlnot4cLRzsKrBw51tnWoyPFVd2PkhLVKhONLHAYB/YDeQGdgP0J26QIR2QksBWYC3wMTVXV6cL0kxFgdeB842le0HHgWWAC0Bi4HDgAOBcaIyBGqmqhfIX4B/lxBnWUJeldylBQ7TZxu+/WM7t7mBzprXfgXBty5zlkHolHbhIZoUkzV2Vl4zK2w1zVwvHq90LEExlSk60nOAqJ+8z6FA89MWzgxKS2Fn9+ANT9Dk47OEglNOsGYoKUyWveFg84N/wy/atWdRfY+vD789fK6p4LtdzCcMxI2/tnpJpv/OexYGzqL66gbyx/Xk+ESnuCISAvgAuBc4DC8rUTlpZd1gQOB7sD5vmdtAD4BXgPGqybl1/vfEkhufgCOV9XNZQGLPAG8B5zgi+1O4OYEvXuDqr6XoGelx7rZ3mbW+q2i32Mkv5oz6n/ppEDZimmW4FRluzc7/wAHz8AAOODUyLMwjImky4neBGf+Z86qxuG6cjLN9JG+jYdd8qoFJRICJ/0d8qJoTel5kTP2ZkuY33vjSUSadITT/uUcq0LRdmdbhl0bnS6zBq1jf2YGSVgXlYgMFJF3cVpAHsFptcnHm9SUl6Bo0HUBmgGXAmOBZSLyJxFJ2A59vnE1d7jef6k7uQFQ1T2+GPy/iv5BRJokKoYqL7h7Ktbf0EPG4eTAaqXZastyeHpAhOTmNGdqsDGxKjzcu+v0ro1V59+Jqc+GlgW3khx6SfStL/kFMOCW8NfKmyIeDRFn7GSjdk6yVMWTG0hAgiMiJ4vINOAL4HQCrULuxEaAImAh8BVOl9DrwCjgbWAM8BOw3lc3uKVHcMbF3AcsF5FHRCQRW5Eeh5NEgdNC9Eu4Sqq6zhcvQA2gnPlyOWb1DO95tN1TfuHG4Ziq6aP/g61Bv1kW1IbTH4dzXy1/fIExkeQXOFu8uM2rAtPF1/0K6+eUX6dGAzjuL7E995DzvUtsgDM2p16L2J6TA+JOcESkr4hMBj7EGZ8SnJjMAP4FnAN0BeqoahdVPVZVz1LVC1X1UlU9R1VPUdVDVbUl0AgYANyIkwhtcb8WqAXcACwSkbtFpDJzToe4jiv6G+O+Xs4w9BwTPE1x/56x3R+8Jsqamd71GEzVMH+cs3CY2349ncGMh16aEauamioseDZVVVgP55d3ved1W0Cd5oFzyYdT/xl9l75fuFac4PWCDBDnGBzf1OWLCU1qvsFp6XhXVVeFubVCqroVmOT7PCYi+TgJz1CcsT2Nfe+sgzMeZoSIXKyqkyI8sjw9XMcVDWh2z2HuEbFWbLqKyCScQcz1gM3AXGA88Iyqrinv5rQrKYK1QY1esbbg1G0ODdvClqXOeWmJsy+VbcZYdewrgbG3e8sKj4BLP3AGRhpTWf5VjdW3o7h/VePgloxMoRqa4Bx3p7N689YVsHG+8+9ek47h76/IIec7a9X8/AY0KHT2iDIh4m3BuZRAcrMdeBQ4QFWPVtUn401uwlHVfar6par+AdgfuBDw7xImQBucrqZ4dHEdL6mg7grA97eLziIJ+ZW0JdAfJ2krAJrjDHi+G1giIhGGy0dPRMJ+EmLdbGfapl/91rH/NgKhrThVpX89lxTtgLeugEcOgE9uhuJdgWvTR8KGua7KAic+aMmNSZzajUN/6Rlza+a29q6bDRvmBc7zqkG3U5yWzIZtoONx8Sc34AywHvYc3LoErvvBmU5vQlRmDM5m4HagUFVvUtW5Fd1QWaparKqvq+pROK06Yyl/ZlZFGrqON1Tw7hLAPz28Gk4LUmUsBP6DM4vrXGA48BCwyHe9Bk4L1l8r+Z7kqWz3lJ+Nw8l8P7wCs95ylpWf+iyMPMlZKn73ZvjyPm/dnhfF//+CMZEE7y4+71N4daiz0W+mCW696Xick6QlWq1GkTfXNHFPE78feCiB68HEzNcldZKIDACaxvkY95rSe6KovxtnjBA4XUo74nzvQFWdGO6CiNyO04LjXx/nLyLyeZxdcMlV2QHGfraiceZb/p33fPUMePZY589ut2viYUEdGHRnSkMzOaL3CJj2PGxeHChbNhleOgUufjtj9j9CFWa94y078Kz0xJLj4mrBUdU/pzO5cVPViar6TsU1M0ek5MZ3bZ+q3omzerLfHZHqp1Vlp4j7tTjI2YTRb/sqp5/aZI5Ni0LLdqyBOR95y47+v8z5QWOyS426cNkYaN7dW752Frx4AmxaHP6+VFvzM2xaGDjPKwhtfTIpket7UblbYKKZjVXLdbw9wbEEu5vAukDHikitcupGpKphP5VWUgRrg1Ywjrdbolr10NYfG4eTOVSj++HRoA0ccW3y4zG5q/5+cNknoRv1bl4Cb47IjK1egrunOg2CWg3TEkquS1qCIyL3J+vZCbTFdVzu4n2+RQH9O0iWEFj4Lyl8A7X9o9RqAO2S+b6YhRtgXCfenkKgTdBA4+WW4GSMnRug2J3PRxj2dvzdUBBXHm5M9Go1gkvehc4neMtXz0j/+L1ws6eseyptktmC8ycR+U8Sn58IrmHuFSYQrXFWZgaYn6RtI4JtdB03ilgrHRI1wNjPZlJlrs1BrTctDoQzn4Z81yypwiOhx7DUxmVyV/XacP5r0OEYb/nM18NWp2g77NqU9LBY/ZP370t+dWcvLZMWye6iulpE/utr/chEs1zHfSLWCr0+K2KtxHK3Km1J0TujEzL+pmflntc6aCbV6hmwbAos/BLmfOK8LxOan3NR8Pibxu2h5wVwxedw8Hlw2FVwwX9tMT+TWvkF0Pcab9msd5zuc7dpL8ADbeCRbs4msMkU0j11vHebCZNSqUg8zgMaiMgw375OcRMRSXDLyVjgJt/xCTh7aEXiXr046ctoish+BNbpKcbZaT1zhMygquQu0fX3c7q5tvkGF+8rdgYOujXtAn2ucBa5Kq9Pu2g7fHiD0wrU62IYGGHvFhOdkASng/Pf/XvC0DB77RiTKp2Oh9pNnP2pAPZsgXljofvpzvnmJTDmFkBhX5GzIOXB5yVnarV1T2WcZLbgTCbQWX8i8JmI1C+nfrlE5FRgZiICc/kSZ/8rgONF5MAI726Ob4dznOnkYXYTTLi7CXz/vlLVpI75iUkiBxi7BY/DCbZhHnx6K/zzAPjguvAze8BJbma95ayO/OV9sGBc5WPLZZESHGPSrVr10K7RmaMDxxMf9m5uuWdraOtzoswbG1iRHSC/hrMTukmbZCY4xwOfEPgh3R+YEOsmmSLSW0S+wEkquldUPxa+xfv8q5QJ8IqIeMa6+Pa6epnAwn5PqKp7bIy77ksior7P3RHqPCAihZFiEpF8EbkHuNpVfF+k+mmx9hfvAOMGbSo3wNiv0+Do6u3dBT+87OxcvTJoh425Y5zkxm36y5WPLZdZgmMy2SHne8/njXXG22xcCDP+F1p/8VeJj0EVJjzgLet2srM7t0mbpHVRqeoeETkDJzm4EGfKc0/gaxEZrKrLy7tfRNriLCh4HoFtIZIxCOM/wDCcLRIOBX4SkWeABTgDi6/A2SsKYDZwbyXf91vgVt9Gpd/gDHTeijMFvTvOqsbunyD3qWoS/kZWQkj31CGJee7B5zp7zMz7DCTPGUhYUNsZ27F0svc3MXBm9rx+EVz1pdPFtWcrfHRj6HPnjnH+wUvGSqK5IDjBadQ+PXEYE87+h0KTzs7+TuD88jXrbWfBUP/eVW6LJ8KAmxMbw7yxof8uHn1T2KomdZI6BkdV9wEXi8hG4A84CUoXYJKIDAm3vYOINMRZxff3QHWSl9j4Yyz2JWJv4exp1YbwScwPwFm+zUArS3BatPqXU2cXcKuqPpGA9yVWomdQ+eUXwJB7nU+w7Wud7QKmj4RtK13lq+H1C521MT6701kkMJj/H7y+VyUmzlyye7N3peL8GlC/VfriMSaYiNOK88XfAmXfPgFbloWvv+w72LsHCspZ+mzvHvjo/2DZt9DyIKcbrMsJ4ZdBUIWvHvSWdTvVuc+kVUoW+lPV6wmMKVGcJOJrEentryMiBSJyE84eTf+Hs/YLvvqCs5v3kCTFtxmnS+184GNgFc7A3rXAFzjdRf1UNcLfmJicgDOw+S2cMUWrfe/ajbOh5ye+620yMrmBxA8wjka9FjDwZrh+Jhz+O++1VT/AK2c43VaRzHgtufFlq+AF/hq1g7xcXx/UZJyDz/Web14CWhq+7r6iipeh+O5p+Om/zpTvXz+AN4fDw53hnWtg4RfeGZ3zPwsd12O7e2eElE3fVtV7RGQ98LivqCnwhYicCeyHM86kkEAS5E9s5gB3qurbSY5PgdG+T7zPGAGMqKDOd8B35dXJaMkaYByt/GpOC8/GBc4/LH7BeyU1LIQtyylr/Fv1I6z7FZofgImBjb8xVUHDQmh3NCz5Ovz1Bm1gq2tUxOKJ0P7oyM+bHWYeSfF2Z52dma87g4dP+zfUbQ4TwrTe7Hdw7F+DSbiU/iqmqv8BLsJZCVhxNqwcB7wKtCWQ3AhOa8aVQI9kJzcmBpuXehd4S9QA41jk5cOwF6BZtwgVBIY+Bx0Geotn/DfpoWWd4BYcS3BMpjr4vPDlLXrAsbd7yxZH3A7QWbm7oplW8z6Fp/rBZ392WpDdbFmKjJHytmZVHQ38Dm8y4x5nswn4I9BFVV9UjdTOaNKiWRe4bTn8fiqc9Uz6/jLXrA8X/M9Ztj1Y36ug8HDoeZG3fOYbsK8ktL6JLNwif8Zkou5neDft9Tv2dmg/wFu28nso2hFaF2DBeDzDPus0g7otQuvt3uyM9XHrekriJl2YSktpgiMijUXkEZxuKv//Qe4BxG8CHVT1n6paFPIAkxny8qFZV2dg36GXpi+Oxh3g3Fcgz9XT2qAQBt3lHHc7FarXC1zbsQYWfZnaGKs666IyVUXN+tDtFG/Zfj2dnbwbtIbGHQPlpSXOSunhLPjce37ocLjxV7j0A2jTr/wYjrGxN5kkJQmOiNQRkb8Ai4AbCOzc7W69EeAkoG+4ZxgTVvsBcMHr0OIgZz+kS96BGnWda9Vrw4FneutbN1VsLMExVckRvwfxbRko+TDkb4EtRIJbcZaE6aYq3edrwXHpPNj5pa7DQLhsjLOprLub3q/rydZ6k2GSmuD4ZkZdj5PY3EVgN25wEppJwO3APpwkpy7wkYgMTWZcJst0Hgy/nQSXj4Gmnb3Xel7oPZ/zsXfas4msaDvsXBc4z6vmjLkyJlO16g0jPnbWuRn+oTepCR5UHG4czqofYbdrU86aDaCVaxvCvHw46v/g6gneaeD5NeCY2xLyJZjESdosKhG5DCepaUNgNWN/S8084E+q+p6v7iyc2Us1caaHjxaR36nqc8mKz+SIwiOcqc2blzjn+4qc/WL6XJ7OqKqG4AHGDQudWWzGZLK2RzifYO2CEpzVPzm/7LjH8c0P6p7qeFz4/+dbHAhXfgE/vuo856BzbOZUBkpmC84LBJIb/zibDcC1wIH+5AZAVT/C6Z7a7qubDzwtIn9KYnwmF4jAIUGtON/8y2mdMOXbbDOoTBap2xyau3b70VJnhXS34PE35W0fU606HHYFnP7v8qecm7RJ1SDjPTjbLnRS1ad8Kxx7qOpEnJWENxBo6blPRB5OUYwmWx1yPoFGRHw7DNtgwArZ+BuTbYLH4Sx2rZuzcwOsDJry3en45MdkkibZCY4CL+FM+f6zqpb7a7Oq/gAMwFkDx5/k3CgiL4qILZ9q4tOoLfS92ls24zWY9U564qkqLMEx2Sa4m8o9DmfhF3gm9bY82FlB3VRZyUwaPgN6qerlqrqywto+vv2pjsIZpwNOkjMceFtEwgxdNyYKg/8aujDghzf4Vjs2YdkifybbtOuPpzV33S+wY71zHDz+pnM53VOmSkhagqOqJ6rqz3HeuxwnyXG3F54OfJqI2EwOKqjlrH7snt5ZtBXevcaZGmpCWQuOyTa1GoVO5f7gWmew8cKg6eHljb8xVULGdvuo6kbgGOArAuvlDCzvHmPK1bIHHP9Xb9nSb2DSo+mJJ5Pt3e3dtV3ynFlUxlR1HY/1ns/7FJ7sB7s2BspqNoDWh6U2LpNwGZvgAKjqDuBE4MN0x2KyRL/fOFM/3SY8AFtXpCeeTOWfVu/XoDVUq5GWUIxJqCOvg7otvWU71nrPOxxrSyJkgYxOcAB8WzYMBUalOxaTBfLy4Mz/QO0mgbLSElgwLn0xZSLrnjLZqnZjuHIc7N8rch0bf5MVMj7BAVDVfap6Kc4eVsZUTr2WcNhV3rI1s9ITS6YKHmDcyDbZNFmkYRu47FPodXH46zY9PCtUiQTHT1VvSHcMJku4l1kHWGsJjoe14JhsV1ATTn8CTvkn5BUEyjsOcn4JMlWedTKa3BSS4PwCqoGN+XKdJTgmF4g4qxHv3wu+exqq13H2sTJZIa4ER0RaqOraimumhog0V9V1Fdc0xqdhIdRo4EwVByjaBluWOvtWGUtwTG5pdSgMfTbdUZgEi7eLaoGI3CMi9RIaTYxEZKCIfAv8Jp1xmCpIxNkwz83G4TiKd8HWoAUQLfEzxlQx8SY4dYA7gGUicr+IpLTDUkROFZGvgS+Avql8t8kiLXt4z20cjmPheGcjQr9G7aB67bSFY4wx8Yg3wSnCWXivAXArsFRE3hSRE5O1Z5SItBKR20VkAfA+cCSBBQB3J+OdJsu1CEpw1sS18Hb2mfOx97zLSemJwxhjKiHeQcbdgEeBM3ESjAKctWqGAptE5ENgDDAx3rE6IpIP9MbZYfwM4DACCY3fOuB2VX0xvi/D5DRrwQm1by/MHeMt63ZKemIxxphKiCvBUdWlwFARORa4FziCQOLRBGdzzOEAIrIImAn8AiwFVgGbgT1ACVATp8urJdAK6AL0AA72lfu5E5vtwGPAI6q6LZ6vwRiad3e2IPB3x2xeAnu2Qc36lXvu5iXw85vQtCsccFp8M7NU4cdRMHM0tOkLx/7ZWaQw2ZZOhj1bAue1GkPhEcl/rzHGJFilpomr6pdAfxE5AbgZp7XFz/+vekegA05rTyz896vreB3wH+BxVd0UT8zGlCmoBU06wYZ5gbJ1s6Hw8PifuWE+PD8I9vhmZ532L+g9IvbnfHEvfP0P53jJ11CzIfS/Lv64ojXnI+9515NsyXpjTJWUkF8JVXWsqh6P0/Lyb5xEJByp4BP28cA44FKgUFX/asmNSZhEjsMp3gmjLwkkNwDTXoj9OZMeDSQ3flOfg9LS8PUTRTV0/I11TxljqqiEtnmr6mzfasOtcFpzHgJ+BPYROYFx8yc664C3gCuB1qo6RFVHqWpxIuM1JmHjcFTho/+D9b96y9fMhB0xLNE09TkYd3do+dZlsGRifLFFa9WP3h3EC2qHbkxqjDFVRFLanlW1FJjg+/xJRGrgjKnpDLQFmgG1fe/fhTOmZjmwBPhZVVeGPNSYZGgRtKJxvGvhTB/pjJcJZ+GXcMh5FT9jxv/gkz9Gvv7jKOhwTFzhRSW4e6rjcU43njHGVEEp6Vz37Qg+zfcxJnMEb9mwbjaU7oO8/OifsepHGHNr5OsLx1ec4MwfB+//zlvmHgANMPsDOHkz1GoUfWyxCO6eOuC05LzHGGNSoEpttmlMwtVrCbWbBM737grdpqA8OzfAG5fCPlfvafBSUAu/KH/8jCp8eqs3mckrgAted7aU8NtXBLPejj62WGxYAOvnBM4lHzoPSc67jDEmBSzBMblNJP6Bxgu/gKePgi3LvOVn/gequ3Yx2bke1pbzzOXfwcYFrpjyYNjz0OUE6Hmxt+6Po6KLLVbB3VPtjoLajZPzLmOMSQFLcIwJ2VncNw6naDt8cjO8fBp8ehvM/9zZp2nvHuf81bNg+2rvvX2vhkPOhw4DveULxkd+f3DScsBpcOCZznHPC/CMz1/1Y3L2zApOcLqdmvh3GGNMCtkCF8aEtODMgn0lMPpiWDTBKVs8EaY8BfnVnS6t4MQGoN3RMORe57jjcd6kYcF4OPrG0HuKd8Iv73nL3K02DQudgcWLvgyUzXgNTnwgyi8uCttWwYqg4XHdTk7c840xJg0yvgVHHF1E5GwRuSfd8ZgsFG6q+Od3BpIbt33FocmN5MHRN8HF70C1Gk5Zp0HeOsunOC1CwX79EIpd5XVbhE7N7hXUTTVzNJQkaMWEOR/D84O9Zfv3ggatE/N8Y4xJk4xqwRGRRjjTyQ/x/fdgoDvgnqv6lzSEZrJZ067OoN7Svc75tpVOa000GhbCWc9C26DtDBq1g8YdYdNC57y0BBZ/HdoyMuM17/nB54WuHNztVKjZILCA4K6NMO9T6H56dDGGs2UZfHILzBsTes26p4wxWSAtLTgiki8iB4rIBSLygIh8LCLLgQ3AFzgbeV4G9MFZL6e8lY6NqZxq1aFZ18jX67eC3pdBw7be8kMugN9MCk1u/IJbcRYGjcPZsszp+nLreVHocwpqwkHnest+fDVyvBWZ/QE82S98ctO4I/S5PP5nG2NMhkh6C46INCfQGuP/HABUD64adK6+T6TrxiROix7hVzGuVhPO/y/s39M537QIVv7gtNC07lP+MzsOgqnPBs6DBxrP+J/3vFVvaN4t/LN6XQzTngucz/8M5nwS+1iZ0lL4+EZnOryb5MMRv4eBt0KNurE90xhjMlDCEhwRKcDpTvInMf5upmbBVYPO/ZtpatB1d73VwFzfx5jEa9nD2fM+2BlPBpIbgMYdnE802h3l7fravNhJkBp3cBKN4O6pcK03fvsd4nxW/xQoe++38JuvvWvlVGTjAmfaulubfnDKP0PHIhljTBUWd4IjIifiHSvTJczzymuVcScy7iTnQ2AGgYRmrqruiDdOY6ISPJMKoP8NcNDZ8T+zRl2n+8rdDbVgPPTtAMsmw5algfL8GtBjaORnicDJj8DIE53xPAB7tsBbV8Bln0B+QXQxrZ7hPS88EkZ8DHkZP9/AGGNiUpkWnE8ITVb81PVf9zX/8QbgF9fnSVedP6hq0MppxiRZ4eHODKYda53zzkNgUALGs3cc5E1wZo52tlr4+U1vvQNOrXgLhjaHwaC7nBlefiumwhd/g8FRTjBcNcN7XtjPkhtjTFZKRBdVuLEx/rItOAnMLFwJjap6tlcWkScxJp0KasHwD+G7Z6BBKzji2tj2o4qk0yAYd1fgfMW00DVnAHpeGN3zjrgWlnztjMHx++Zfzho8nQdHvs8vuAVnv57RvdcYY6oYUdWKa4W7UaQUb9fSTOB7XAmNqoZZDS3is/A9p7214FSOiJT9ocb752sSRBUe6RpoGQqn3v7wf7OiT6h2bnS2iNi+KlBWuwn85huov1/k+0pL4cFC77o71//kDJg2xpgqSCTQxqKqngaXRLRN+5OcpsBk4F+qOi7a5MaYrCYCh/+2/Dp9r4qttahOEzj7BWfmk9+ujfDlfeXft2mhN7mp2TB06rsxxmSJyrTg/AO4Fme6t/shs4BbVHVsDM+yFpwEshacDFO6D+Z+Aku+cRKMou1QtMNZFbnd0TDgj/F1h038hzP+xk/y4dpp0KRj+Poz34R3rgycdzgGLn0/9vcaY0yGKK8FJ+4Ex/fg9sDDQPD0DwXGAzer6k8hN4Y+xxKcBLIEJ0fsK4Gn+nl3Ij/4fBj6TPj6Y++Ab58InPe/AQb/NakhGmNMMiWti0pVF6vq2cAAnPE3EOiyGgRMF5GXRaRNZd5jjAkjvxoM/JO3bOZoWB9huajgGVTu9X2MMSbLJGR+qKpOUtW+wKXAiqDnXwzM9W3JUD8R7zPG+PQYCs0OcBUoTAiz03hpqXeRQLAZVMaYrJbQBTBUdRTQFbgT2Om6VBO4BVgoIteJSEZt8mlMlZWXD8fe5i375V1YE7TtxKZFoQOMbfaUMSaLJXyFL1Xdo6r3AZ2AF4BS1+UmOBtpzhGRSiwRa4wp0+00aHmQtyy4FSdk/ZtDnBlexhiTpZK2hKmqrlPVq4BDgXF4t2ToAIwWkckiclSyYjAmJ+TlwbF3eMvmfASrfgycu4/Bxt8YY7Je0tdoV9WfVXUIcBowx1+Mk+j0A77CO808LcRxnoh8JCIrRKRIRFaLyHgRuTIZ3WoiUkdEbhSRb0RknYjsEZGlIjJaRE5I9PtMFutyorMbudsXrnVxbPyNMSbHpGwTGlX9GGdTzj8AG12XgtvJo9w1MHFEpBFOK9PrwClAK5z1fVoCxwHPAd+JSAzbNlf4zl44qz8/AhyJs+t6DaAQOBf4VEReE5HqiXqnyWIicOzt3rIFnzube4YbYGwtOMaYLJfSXfZUdZ+qPokzPucfQLH7Mk6yM0lErhWRlCQ6vgTifZxEBmA5ziDpC4CbgV995YcCYxIxE0xE2gJjcLrqAKYC1wEXAg8QSAAvxBnHZEzFOg5ydgd3G3s7bJwPRdsCZTUbQKP2qY3NGGNSrFIL/VX65SLtgIeAs/HuTK44icbdwMuaxCBF5HrgMd/pD8DxqrrZdb0m8B7g7zL6h6reXMl3vguc6Tt9EbhKVUtd1wuBr3FacwBO9bWARft8W+gvV638AZ471lvW9ihYOilw3n6As7GoMcZUcUlbyThRRORI4J9AXwItOfiO5wF3qupbSXhvNWAVTveQAgep6i9h6jUHFgF1gCKglapuDK4X5TsPAWb4TpcBXVV1T5h6JwP+pOZ7VT0shndYgpPL3vsdzHgt8vUjr4Mhf4t83Rhjqohkb7ZZaao6WVUPx1kU0L9Ngz/R6Yoz4+p7ETkxwa8+Die5ARgfLrnxxbcOZ3wOOONkzqjEO89zHT8bLrnxGQP41+DvIyIRNhgyJsigv0D1upGv2/gbY0wOyIgEx09V/wt0A+4AdviLcRKdQ4GPEvzKIa7jTyuo675emUQrqnf6uuXcG5barCoTnXot4eibIl+3GVTGmByQcSsKq2oR8ICIPA/cC1yOk4gJoTOuKquH63h6BXW/dx33iFirHCKSB/jX1S8BKtqINBHvDFseS9dVuz9FPfzHZIgatOfz6s0ozFvvKd+mtTn44V8JrNhgjDHJt+TBU1L+zoxqwXFT1fWqeg3QC/g8Sa/p4jpeUkHdFcA+33FniZQ5lK81UNt3vFJVSyqov9R13CViLWOCFFGd+0ouCimfVdqOxP+eYIwxmSdjExw/VZ2lqifirE8zO8GPb+g63lBBHCWAf65tNZwBx0l7n497IHPDSJWMCWds6WFM3tfdUzZTbSiXMSY3ZHyC46eqY4BDEvxY90jMSIN93Xa7jutVgfeZnCbcWXIZm9T53263VmfUvkFpjskYY1Ij48bglMe9VkwWSMn87URME09H36lJoK0nw/Kp1Gp7JJPqtUx3NMYYkxJVKsFJgh1AI99xTQIztyKp5TreHuf7wj0rWe8zBhq0dj7GGJNDqkwXVZJscR03Ka+ib1FA/zYNJcDOZL4vTJ0tkSoZY4wxxivXE5x5ruN2FdRtDeT7jufHuX3ECmCX/3lR7FDe1nU8L2ItY4wxxnjkeoIzy3Xcp4K67uuzItYqh28MkX/zzmpUPGi60u80xhhjclGuJzixrBTsXr24olWPK/1O3zo77utjI9U1VZ+IeD6marE/v6rN/vyyU0ZstpkucW62uQdobZttmkQK/kfV/tyqFvvzq9rsz6/qyvjNNtPFt3jffb5TAV4RkUbuOiJSE3iZwMJ+T0RKbkTkJRFR3+fuCO/8CXjPd1oIPOHbwsH9nELgP66isM8yxhhjTHi5Pk0cnERiGHA0zoaeP4nIMzg7ebcGriCwf9RsnP2xKusG4Aighe/5PUTkVZyViw8CriEwg+o1VbXNoIwxxpgY5HyCo6rFInIG8BZwHNCG8EnMD8BZqro1Ae9cKiIn+d7ZAejn+wT7L85mo8YYY4yJQU53Ufmp6mbgeOB8nHEvq4BiYC3wBXA10E9VlyXwnT8CBwM3Ad/i7E1VBCwH3gROUtWLVLU4Ue80xhhjckVODzLOVjbIuOqxQY5Vm/35VW3251d1lTfI2BKcLOROcIwxxphcYLOojDHGGJP1LMExxhhjTNaxLipjjDHGZB1rwTHGGGNM1rEExxhjjDFZxxIcY4wxxmQdS3CMMcYYk3UswTHGGGNM1rEExxhjjDFZxxIcY4wxxmQdS3CMMcYYk3UswTHGGGNM1rEEx5g0EpFaInKaiDwqIpNEZJ2IFIvINhH5VURGisigKJ91t4hoDJ92Sf7yspo4zhORj0RkhYgUichqERkvIleKSLV0x5hrRKSBiJwrIv8Rke9EZKOI7BWRzSLyk4g8JSKHRfGcl2L5u5SKr83EzrZqMCZNROQi4GmgbhTVPwUuVdX15TzvbuCuGEJor6pLYqhvfESkEfAWcFw51X4AzlLVZamJKreJyC3APUCNKKqPAq5R1V0RnvUSMDzadwfvYm0yg/2GYUz6tCeQ3KwGPgemAeuAOsDRwAVATeBEYJyIHBHpH+UgdwKzKqizLp6gc52IVAfex/nzAVgOPAssAFoDlwMHAIcCY3x/ZtvSEWuO6UIguVkEjANmABuARsAgYBiQD1wMNBeRk1S1tILnXoP9XamSLMExJr2+AR4ExqjqvqBrI0XkHzj/UO8HHAzcSnStNJNUdUIiAzVlfksgufkBOF5VN/svisgTwHvACUB3nGTz5hTHmIsU+Bh4GJiood0Tz4rI0cAnOL9YDMFppRlZwXM/s5bOqsnG4BiTPk+q6lGq+lGY5AYAVZ0NXO0qGpGSyExYvnE1d/hOFafbcLO7jqruAS4FdvqK/iAiTVIXZc66RVVPVdWvwiQ3AKjq18BtrqIRKYnMpIUlOMakSfAPxnKMIfDDslBE6icpJFOx44BmvuPxqvpLuEqqug543XdaAzgjBbHltBj+Pr3pOj4oGbGYzGAJjjEZzte64x53UytdsRiGuI4/raCu+/qJSYjFxGe769j+LmUxS3CMyXAi0pxAq8EuIOJMKpd7RGShiOzxTTmfLyKv+qak24yP+PVwHU+voO73Ee4z6eX+s1gaRf3nRGSZbxmALSIyW0SeE5EByQrQJIYNMjYm87nH4HwaxawPCAyCBaeLpB7QCWf2yHcicp6qRvOPu/Hq4jpeUkHdFcA+nFk7nUVEIo0NMSnl/vv0cRT1j3cdVwca4MySu1JEPsYZh7UpgfGZBLEEx5gMJiId8A6KfLCCW/YAXwLfAouBYqAlMBBnHEg+0A+YLCJ9VXVlwoPObg1dxxvKq6iqJSKyDWeKcjWcqf87kheaqYiIHAlc5jvdAzxWTvXtOEs3TMVZCmAfzjIAQwh0VZ4CfCUi/W0pgMxjC/0Zk6FEpA4wGWd6ODizrq4tp34PYFWk3yZF5ECc6cudfEXjVHVw4iLOfiJSDBT4TgtUtaSC+iuB/X2n+6vq6mTGZyITkZY43YatfEU3q+o/ItTtDcxR1Z0Rrh+Ns9Bjc1/R86p6VYJDNpVkCY4x5RCRK3F+a6s0Vb07hvfmA+8Ap/uKfgD6+6Ygx01EOgI/ExhcebiqfleZZ+YSS3CqJt8vC18AfX1FHwOnVabLUET6A18DgtO609ZaRDOLJTjGlENEpuB06VRatMu5i0ge8DLOeBmAucAA39TjShORp3FWZwX4m6r+JRHPzQUisgmnywmgnqqW2+UUa32TeCJSEyeh8W+r8Q1wQqTWmRif/SnOgo4AV6jqi5V9pkkcm0VlTAbxzXB6hkBysxAYlKjkxmeC67hrAp+bC7a4jstdvM+3KKB/zaISAmsZmRTxbavxDoHkZipwciKSG58JrmP7u5RhbJCxMeVQ1cNT/MongCt9x0uB45LQ7L3RddwoYi0TzjycPcQA2lH+NOPWOIO6AebbDKrUEpECnEX9TvIV/QicmODBwPZ3KYNZC44xGUJEHgN+5ztdgZPcJGMnanfLw5YkPD+buTcw7VNBXff1ijY+NQnkaz37H4ExbD8Dg2NY7Tha9ncpg1mCY0wGEJGHgOt9p6txkptFSXrdQNfxvCS9I1uNdR2fELGWw716cUWrHpsE8Q3QfxVn53CA2Tgbom6MfFfc7O9SBrNBxsakmYjcS2ADx7XAMao6J0nvCp5F1V9VJyfjXdnI1zKwCmdlaQUOCrcflW/16UU4a9/sAVon6QescfEN0B+Js9kpOEnHQFVdk4R3HQlMwplFVQq0U9XliX6PiZ+14BiTRiLyZwLJzXqcAcUxJze+LRjO9v32GqnOgTgtCf7kZoIlN7HxTQu/z3cqwCsi4hl74Zu18zJOcgPwhCU3yecaoO9PbhYAx8aa3IjIpSIyuLwtTUTkKOBdnP8HAF6x5CbzWAuOMWkiIlcBz7qK7gJmRnHrJFX1rKIrIjcAj+IkSZ8CM4A1wF6gBXAMzkrG/okFq4AjkjTGJ6v5ZuaMI7AdxnKcH6wLcAYWX4GzlD843SNHqurWVMeZa0TkfgKrfu8FbsL5s6nIZ6patpmtbyzc9b57x+K0eK4ndCVjf3LzC3CUqm6p9BdhEsoSHGPSREReAobHceuxqjoh6Fk34CQ40fgauMT2ooqfr9XmLQLTj8P5ATjLksjUEJEJeMfERKu9qi5xPecxAuPhKvIucJW10GUmmyZuTHZ4FVgJHIEze2d/oClON8k2nN9GpwCvBydHJnaqullEjgfOBS4BeuF8vzfj/Eb/OjCyopWOTUZ6GGen+MOBQ3FaQJsCNYGtOHu8TcbplvohXUGailkLjjHGGGOyjg0yNsYYY0zWsQTHGGOMMVnHEhxjjDHGZB1LcIwxxhiTdSzBMcYYY0zWsQTHGGOMMVnHEhxjjDHGZB1LcIwxxhiTdSzBMcYYY0zWsQTHGGOMMVnHEhxjjDHGZB1LcIwxxhiTdSzBMcYYY0zWsQTHGGOMMVnHEhxjjDHGZB1LcIwxxhiTdSzBMcYYY0zWsQTHGGOMMVnHEhxjjImBiIwVEXV9rojyvloiMino3geTHa8xuUpUNd0xGGNMlSEiBwM/EvgFcTnQWVWLyrmnGvAecIqr+DlVvTpZcRqT66wFxxhjYqCqM4FXXEVtgN9Fqi8iAryEN7l5A/hNMuIzxjisBccYY2IkIq2BeUAtX9EGoIOqbg9T91/Ada6iscDpqlqc9ECNyWHWgmOMMTFS1RXAo66ipsBNwfVE5C94k5vJwFBLboxJPmvBMcaYOIhIfWAB0MxXtB2nFWeD7/pvgadct8wEBqrqllTGaUyushYcY4yJg6puA/7qKqoH3AEgIucBT7iuLQCGWHJjTOpYC44xxsTJNzvqF6CLr6gIuBan5abAV7YSOEpVl6Q8QGNymCU4xhhTCSJyFvBOhMsbgQGqOjuFIRljsATHGGMqTUS+Bo4KKt4BDFLVqWkIyZicZ2NwjDGm8j4JU3aWJTfGpI8lOMYYUwkicipwT5hLLVMdizEmwLqojDEmTiIyAPiUwIJ/bsuALuVt4WCMSR5rwTHGmDiIyKHAhwSSm53AW64qhXgX+TPGpJC14BhjTIxEpCvwNYFF/oqB04ApwCKgia98C9BRVTelOkZjcp214BhjTAxEpA3wGYHkphS4WFU/8y3+94CrekPgz6mN0BgD1oJjjDFRE5GmOC033VzFv1XVp111auJsxNnGV1QMdFPVxSkL1BhjLTjGGBMNEamHM6DYndz82Z3cAKjqHrxbOFQH7k9+hMYYN2vBMcaYCvhaZcYAx7iKH1XVGyPUzwdmEUiGFOinqtOSGacxJsBacIwxphy+/aZG401uXgFuinSPqu4D7nQ/Bng4GfEZY8KzFhxjjIlARAR4GbjEVfwhMFRVS6K4fxrQx1V0uqp+mNgojTHhWIJjjDERiMi/8K5lMxE4wTfOJpr7jwc+dxXNBg72tfAYY5LIEhxjjDHGZB0bg2OMMcaYrGMJjjHGGGOyjiU4xhhjjMk6luAYY4wxJutYgmOMMcaYrGMJjjHGGGOyjiU4xhhjjMk6luAYY4wxJutYgmOMMcaYrGMJjjHGGGOyjiU4xhhjjMk6luAYY4wxJutYgmOMMcaYrGMJjjHGGGOyzv8DJTwhc4W+LiQAAAAASUVORK5CYII=\n",
-      "text/plain": [
-       "<Figure size 576x432 with 1 Axes>"
-      ]
-     },
-     "metadata": {
-      "needs_background": "light"
-     },
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "fig, axs = plt.subplots(1, figsize=(8,6))\n",
     "\n",
-    "axs.plot(grid, pmf2, linewidth=4, label=r'eABF/MBAR')\n",
+    "#axs.plot(grid, pmf2, linewidth=4, label=r'eABF/MBAR')\n",
     "axs.plot(grid+bin_width, pmf_ti2, linewidth=4, label='eABF/CZAR')\n",
     "\n",
     "axs.tick_params(axis='y',length=8,width=4,labelsize=30, pad=10, direction='in')\n",
diff --git a/setup.py b/setup.py
index 62ae495..806d1e0 100644
--- a/setup.py
+++ b/setup.py
@@ -10,9 +10,10 @@
         "adaptive_sampling",
         "adaptive_sampling.colvars",
         "adaptive_sampling.processing_tools",
-        "adaptive_sampling.sampling_tools"
+        "adaptive_sampling.sampling_tools",
+        "adaptive_sampling.interface",
     ],
-    version="1.0.0",
+    version="2.0.3",
     license="MIT",
     description="Sampling algorithms for molecular transitions",
     long_description=open("README.md").read(),
@@ -20,7 +21,7 @@
     author="Andreas Hulm",
     author_email="andreas.hulm@cup.uni-muenchen.de",
     url="https://github.com/ochsenfeld-lab/adaptive_sampling",
-    download_url="https://github.com/ochsenfeld-lab/adaptive_sampling/archive/refs/tags/v1.0.0.zip",
+    download_url="https://github.com/ochsenfeld-lab/adaptive_sampling/archive/refs/tags/v2.0.3.zip",
     keywords=[
         "computational chemistry",
         "molecular dynamics",
@@ -33,9 +34,9 @@
         "scipy>=1.7.0",
     ],
     setup_requires=["pytest"],
-    tests_requires=["pytest"],
+    #test_requires=["pytest"],
     classifiers=[
-        "Development Status :: 5 - Production/Stable",  # Chose either "3 - Alpha", "4 - Beta" or "5 - Production/Stable" as the current state of your package
+        "Development Status :: 4 - Beta",  # Chose either "3 - Alpha", "4 - Beta" or "5 - Production/Stable" as the current state of your package
         "Environment :: Console",
         "License :: OSI Approved :: MIT License",
         "Programming Language :: Python :: 3.8",


Test patch:
diff --git a/adaptive_sampling/test/colvars/resources/traj.xyz b/adaptive_sampling/test/colvars/resources/traj.xyz
deleted file mode 100644
index fcb0463..0000000
--- a/adaptive_sampling/test/colvars/resources/traj.xyz
+++ /dev/null
@@ -1,31 +0,0 @@
-4
-
-H	0	0	0
-H	1	0	0
-H	0	1	0
-H	0	0	1
-4
-
-H	0	0	0
-H	2	0	0
-H	0	1	0
-H	0	0	1
-4
-
-H	0	0	0
-H	3	0	0
-H	0	1	0
-H	0	0	1
-4
-
-H	0	0	0
-H	4	0	0
-H	0	1	0
-H	0	0	1
-4
-
-H	0	0	0
-H	5	0	0
-H	0	1	0
-H	0	0	1
-
diff --git a/adaptive_sampling/test/colvars/rmsd_test.py b/adaptive_sampling/test/colvars/rmsd_test.py
index 38fb3c1..01ab5d2 100644
--- a/adaptive_sampling/test/colvars/rmsd_test.py
+++ b/adaptive_sampling/test/colvars/rmsd_test.py
@@ -26,85 +26,6 @@ def get_sampling_data(self) -> SamplingData:
             0.0,
         )
 
-
-@pytest.mark.parametrize(
-    "a, b, expected_nofit, expected_fit",
-    [
-        (
-            torch.tensor(
-                [
-                    [0, 0, 0],
-                    [0, 0, 1],
-                    [0, 0, 2],
-                    [9, 9, 7],
-                    [9, 9, 8],
-                    [9, 9, 9],
-                ]
-            ),
-            torch.tensor(
-                [
-                    [30.50347534, -20.16089091, -7.42752623],
-                    [30.77704903, -21.02339348, -7.27823201],
-                    [31.3215374, -21.99452332, -7.15703548],
-                    [42.05988643, -23.50924264, -15.59516355],
-                    [42.27217891, -24.36478643, -15.59064995],
-                    [42.66080502, -25.27318759, -15.386241],
-                ]
-            ),
-            26.6020,
-            0.0979,
-        )
-    ],
-)
-def test_kearsley(a, b, expected_nofit, expected_fit):
-    k = Kearsley()
-    a = torch.flatten(a).float()
-    b = torch.flatten(b).float()
-
-    rmsd_nofit = rmsd(a, b)
-    rmsd_kearsley = k.fit(a, b)
-
-    assert float(rmsd_nofit) == pytest.approx(expected_nofit, rel=1.0e-3)
-    assert float(rmsd_kearsley) == pytest.approx(expected_fit, rel=1.0e-3)
-
-
-@pytest.mark.parametrize(
-    "a, b, expected_nofit, expected_fit",
-    [
-        (
-            torch.tensor(
-                [
-                    [0, 0, 0],
-                    [1, 0, 0],
-                    [0, 1, 0],
-                    [0, 0, 1],
-                ]
-            ),
-            torch.tensor(
-                [
-                    [0, 0, 0],
-                    [0, 0, 1],
-                    [1, 0, 0],
-                    [0, 1, 0],
-                ]
-            ),
-            0.707,
-            0.0000,
-        )
-    ],
-)
-def test_kearsley_1(a, b, expected_nofit, expected_fit):
-    k = Kearsley()
-    a = torch.flatten(a).float()
-    b = torch.flatten(b).float()
-
-    rmsd_nofit = rmsd(a, b)
-    rmsd_kearsley = k.fit(a, b)
-
-    assert float(rmsd_nofit) == pytest.approx(expected_nofit, abs=1.0e-3)
-    assert float(rmsd_kearsley) == pytest.approx(expected_fit, abs=1.0e-3)
-
-
 @pytest.mark.parametrize(
     "a, b, expected_nofit, expected_fit",
     [
@@ -257,33 +178,6 @@ def test_quaternion_1(a, b, expected_nofit, expected_fit):
     assert float(rmsd_quaternion) == pytest.approx(expected_fit, abs=1.0e-3)
 
 
-@pytest.mark.parametrize(
-    "images, interpolated_images",
-    [
-        (
-            [
-                torch.tensor([0.0, 0, 0, 0, 0, 0]),
-                torch.tensor([1.0, 0, 0, 0, 0, 0]),
-            ],
-            torch.tensor(
-                [
-                    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [0.2, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [0.4, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [0.6, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [0.8, 0.0, 0.0, 0.0, 0.0, 0.0],
-                    [1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
-                ],
-            ),
-        )
-    ],
-)
-def test_interpolate(images, interpolated_images):
-    i = interpolate_coordinates(images, n_interpol=5)
-    i = torch.stack(i)
-    assert torch.allclose(i, interpolated_images, rtol=1.0e-5)
-
-
 @pytest.mark.parametrize(
     "file, coords",
     [
@@ -298,30 +192,6 @@ def test_read_xyz(file, coords):
     assert torch.allclose(xyz, coords.float())
 
 
-@pytest.mark.parametrize(
-    "file, coords",
-    [
-        (
-            "resources/traj.xyz",
-            torch.tensor(
-                [
-                    [0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1],
-                    [0, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 1],
-                    [0, 0, 0, 3, 0, 0, 0, 1, 0, 0, 0, 1],
-                    [0, 0, 0, 4, 0, 0, 0, 1, 0, 0, 0, 1],
-                    [0, 0, 0, 5, 0, 0, 0, 1, 0, 0, 0, 1],
-                ]
-            )
-            / BOHR_to_ANGSTROM,
-        )
-    ],
-)
-def test_read_traj(file, coords):
-    xyz = read_traj(file)
-    xyz = torch.stack(xyz)
-    assert torch.allclose(xyz, coords.float())
-
-
 @pytest.mark.parametrize(
     "coords1, coords2, rmsd_fit", [("resources/m1.xyz", "resources/m2.xyz", 0.000)]
 )
@@ -338,54 +208,3 @@ def test_rmsd(coords1, coords2, rmsd_fit):
     assert np.allclose(grad1, grad2)
     assert np.allclose(grad2, grad3)
 
-
-#@pytest.mark.parametrize(
-#    "coords, internals",
-#    [
-#        (
-#            torch.tensor([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]),
-#            torch.tensor(
-#                [
-#                    1.0000,
-#                    1.4142,
-#                    1.4142,
-#                    1.0000,
-#                    1.4142,
-#                    1.4142,
-#                    1.0000,
-#                    1.4142,
-#                    1.4142,
-#                    1.0000,
-#                    1.0000,
-#                    1.0000,
-#                ]
-#            ),
-#        )
-#    ],
-#)
-#def test_internals(coords, internals):
-#    coords = coords.float()
-#    close_indices = find_closest_points(coords)
-#    ints = get_internal_coords(coords, close_indices)
-#    assert torch.allclose(ints, internals)
-
-
-@pytest.mark.parametrize(
-    "coords, close_points",
-    [
-        (
-            torch.tensor([0, 0, 0, 1, 0, 0, 2, 0, 0, 3, 0, 0]),
-            [
-                np.array([1, 2, 3]),
-                np.array([2, 0, 3]),
-                np.array([1, 3, 0]),
-                np.array([2, 1, 0]),
-            ],
-        )
-    ],
-)
-def test_find_close_points(coords, close_points):
-    points = find_closest_points(coords)
-    close_points = np.stack(close_points)
-    points = np.stack(points)
-    assert np.allclose(points, close_points)