def interpolate(self):
        """
        Finds a set of n_images from self.s1 to self.s2, where all sites except
        for the ones given in relax_sites, the interpolation is linear (as in
        pymatgen.core.structure.interpolate), and for the site indices given
        in relax_sites, the path is relaxed by the elastic band method within
        the static potential V.
        """
        images = self.__s1.interpolate(self.__s2, nimages=self.__n_images,
                                       interpolate_lattices=False)
        for site_i in self.__relax_sites:
            start_f = images[0].sites[site_i].frac_coords
            end_f = images[-1].sites[site_i].frac_coords

            path = NEBPathfinder.string_relax(
                NEBPathfinder.__f2d(start_f, self.__v),
                NEBPathfinder.__f2d(end_f, self.__v),
                self.__v, n_images=(self.__n_images + 1),
                dr=[self.__s1.lattice.a / self.__v.shape[0],
                    self.__s1.lattice.b / self.__v.shape[1],
                    self.__s1.lattice.c / self.__v.shape[2]])
            for image_i, image in enumerate(images):
                image.translate_sites(site_i,
                                      NEBPathfinder.__d2f(path[image_i],
                                                          self.__v) -
                                      image.sites[site_i].frac_coords,
                                      frac_coords=True, to_unit_cell=True)
        self.__images = images