Implements path relaxation via the elastic band method. In general, the
        method is to define a path by a set of points (images) connected with
        bands with some elasticity constant k. The images then relax along the
        forces found in the potential field V, counterbalanced by the elastic
        response of the elastic band. In general the endpoints of the band can
        be allowed to relax also to their local minima, but in this calculation
        they are kept fixed.

        Args:
            start, end: Endpoints of the path calculation given in discrete
                coordinates with respect to the grid in V
            V: potential field through which to calculate the path
            n_images: number of images used to define the path. In general
                anywhere from 20 to 40 seems to be good.
            dr: Conversion ratio from discrete coordinates to real coordinates
                for each of the three coordinate vectors
            h: Step size for the relaxation. h = 0.1 works reliably, but is
                slow. h=10 diverges with large gradients but for the types of
                gradients seen in CHGCARs, works pretty reliably
            k: Elastic constant for the band (in real units, not discrete)
            min_iter, max_iter: Number of optimization steps the string will
                take before exiting (even if unconverged)
            max_tol: Convergence threshold such that if the string moves by
                less than max_tol in a step, and at least min_iter steps have
                passed, the algorithm will terminate. Depends strongly on the
                size of the gradients in V, but 5e-6 works reasonably well for
                CHGCARs.

def string_relax(start, end, V, n_images=25, dr=None, h=3.0, k=0.17,
                     min_iter=100, max_iter=10000, max_tol=5e-6):
        """
        Implements path relaxation via the elastic band method. In general, the
        method is to define a path by a set of points (images) connected with
        bands with some elasticity constant k. The images then relax along the
        forces found in the potential field V, counterbalanced by the elastic
        response of the elastic band. In general the endpoints of the band can
        be allowed to relax also to their local minima, but in this calculation
        they are kept fixed.

        Args:
            start, end: Endpoints of the path calculation given in discrete
                coordinates with respect to the grid in V
            V: potential field through which to calculate the path
            n_images: number of images used to define the path. In general
                anywhere from 20 to 40 seems to be good.
            dr: Conversion ratio from discrete coordinates to real coordinates
                for each of the three coordinate vectors
            h: Step size for the relaxation. h = 0.1 works reliably, but is
                slow. h=10 diverges with large gradients but for the types of
                gradients seen in CHGCARs, works pretty reliably
            k: Elastic constant for the band (in real units, not discrete)
            min_iter, max_iter: Number of optimization steps the string will
                take before exiting (even if unconverged)
            max_tol: Convergence threshold such that if the string moves by
                less than max_tol in a step, and at least min_iter steps have
                passed, the algorithm will terminate. Depends strongly on the
                size of the gradients in V, but 5e-6 works reasonably well for
                CHGCARs.
        """
        #
        # This code is based on the MATLAB example provided by
        # Prof. Eric Vanden-Eijnden of NYU
        # (http://www.cims.nyu.edu/~eve2/main.htm)
        #

        # print("Getting path from {} to {} (coords wrt V grid)".format(start, end))
        
        # Set parameters
        if not dr:
            dr = np.array(
                [1.0 / V.shape[0], 1.0 / V.shape[1], 1.0 / V.shape[2]])
        else:
            dr = np.array(dr, dtype=float)
        keff = k * dr * n_images
        h0 = h

        # Initialize string
        g1 = np.linspace(0, 1, n_images)
        s0 = start
        s1 = end
        s = np.array([g * (s1 - s0) for g in g1]) + s0
        ds = s - np.roll(s, 1, axis=0)
        ds[0] = (ds[0] - ds[0])
        ls = np.cumsum(la.norm(ds, axis=1))
        ls = ls / ls[-1]
        fi = interp1d(ls, s, axis=0)
        s = fi(g1)

        # Evaluate initial distances (for elastic equilibrium)
        ds0_plus = s - np.roll(s, 1, axis=0)
        ds0_minus = s - np.roll(s, -1, axis=0)
        ds0_plus[0] = (ds0_plus[0] - ds0_plus[0])
        ds0_minus[-1] = (ds0_minus[-1] - ds0_minus[-1])

        # Evaluate potential gradient outside the loop, as potential does not
        # change per step in this approximation.
        dV = np.gradient(V)

        # Evolve string
        for step in range(0, max_iter):
            if step > min_iter:
                # Gradually decay step size to prevent oscillations
                h = h0 * np.exp(-2.0 * (step - min_iter) / max_iter)
            else:
                h = h0
            # Calculate forces acting on string
            d = V.shape
            s0 = s
            edV = np.array([[dV[0][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
                                 int(pt[2]) % d[2]] / dr[0],
                             dV[1][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
                                 int(pt[2]) % d[2]] / dr[0],
                             dV[2][int(pt[0]) % d[0]][int(pt[1]) % d[1]][
                                 int(pt[2]) % d[2]] / dr[0]] for pt in s])
            # if(step % 100 == 0):
            #    print(edV)

            # Update according to force due to potential and string elasticity
            ds_plus = s - np.roll(s, 1, axis=0)
            ds_minus = s - np.roll(s, -1, axis=0)
            ds_plus[0] = (ds_plus[0] - ds_plus[0])
            ds_minus[-1] = (ds_minus[-1] - ds_minus[-1])
            Fpot = edV
            Fel = keff * (la.norm(ds_plus) - la.norm(ds0_plus)) * (
                ds_plus / la.norm(ds_plus))
            Fel += keff * (la.norm(ds_minus) - la.norm(ds0_minus)) * (
                ds_minus / la.norm(ds_minus))
            s -= h * (Fpot + Fel)

            # Fix endpoints
            s[0] = s0[0]
            s[-1] = s0[-1]

            # Reparametrize string
            ds = s - np.roll(s, 1, axis=0)
            ds[0] = (ds[0] - ds[0])
            ls = np.cumsum(la.norm(ds, axis=1))
            ls = ls / ls[-1]
            fi = interp1d(ls, s, axis=0)
            s = fi(g1)

            tol = la.norm((s - s0) * dr) / n_images / h

            if tol > 1e10:
                raise ValueError(
                    "Pathfinding failed, path diverged! Consider reducing h to "
                    "avoid divergence.")

            if step > min_iter and tol < max_tol:
                print("Converged at step {}".format(step))
                break

            if step % 100 == 0:
                print("Step {} - ds = {}".format(step, tol))
        return s