diff --git a/ase/neb.py b/ase/neb.py
index 940a8c7e1..69b4dbb6e 100644
--- a/ase/neb.py
+++ b/ase/neb.py
@@ -49,12 +49,15 @@ class Spring:
 
 
 class NEBState:
-    def __init__(self, neb, images, energies, precon):
+    def __init__(self, neb, images, energies):
         self.neb = neb
         self.images = images
         self.energies = energies
-        self.precon = precon
-            
+
+    @property
+    def precon(self):
+        return self.neb.precon
+
     def spring(self, i):
         return Spring(self.images[i], self.images[i + 1],
                       self.energies[i], self.energies[i + 1],
@@ -166,7 +169,7 @@ class FullSpringMethod(NEBMethod):
     """
     Elastic band method. The full spring force is included.
     """
-    
+
     def get_tangent(self, state, spring1, spring2, i):
         # Tangents are bisections of spring-directions
         # (formula C8 of paper III)
@@ -215,12 +218,12 @@ class BaseSplineMethod(NEBMethod):
 
     def add_image_force(self, state, tangential_force, tangent, imgforce,
                         spring1, spring2, i):
-                
+
         # update preconditioner and apply to image force
         precon_imgforce, _ = state.precon[i].apply(imgforce.reshape(-1),
                                                    state.images[i])
         imgforce[:] = precon_imgforce.reshape(-1, 3)
-        
+
         # project out tangential component (Eqs 6 and 7 in Paper IV)
         imgforce -= tangential_force * tangent
 
@@ -240,7 +243,7 @@ class SplineMethod(BaseSplineMethod):
                         spring1, spring2, i):
         super().add_image_force(state, tangential_force,
                                 tangent, imgforce, spring1, spring2, i)
-        
+
         reaction_coordinate, _, x, dx, d2x_ds2 = state.spline
 
         # Definition following Eq. 9 in Paper IV
@@ -259,10 +262,14 @@ class StringMethod(BaseSplineMethod):
     def adjust_positions(self):
         # fit cubic spline to positions, reinterpolate to equispace images
         # note this use the preconditioned distance metric.
-        s, _, x, _, _ = self.neb.spline_fit()
+        _, _, x, _, _ =
+
+        fit = self.neb.spline_fit()
+
         new_s = np.linspace(0.0, 1.0, self.neb.nimages)
-        new_positions = x(new_s[1:-1]).reshape(-1, 3)
-        self.neb.set_positions(new_positions)
+        new_positions = fit.x(new_s[1:-1]).reshape(-1, 3)
+        return new_positions
+        #self.neb.set_positions(new_positions)
 
 
 def get_neb_method(neb, method):
@@ -284,7 +291,7 @@ class BaseNEB:
     def __init__(self, images, k=0.1, climb=False, parallel=False,
                  remove_rotation_and_translation=False, world=None,
                  method='aseneb', allow_shared_calculator=False, precon=None):
-        
+
         self.images = images
         self.climb = climb
         self.parallel = parallel
@@ -460,7 +467,7 @@ class BaseNEB:
         self.energies = energies
         self.real_forces = np.zeros((self.nimages, self.natoms, 3))
         self.real_forces[1:-1] = forces
-        
+
         # if this is the first force call, we need to build the preconditioners
         if self.precon is None or isinstance(self.precon, str):
             self.precon = make_precon_images(self.precon, images)
@@ -605,13 +612,13 @@ class BaseNEB:
                 raise ValueError(f'unknown norm {norm} in spline_fit()')
 
         s = d_P.cumsum() / d_P.sum()  # Eq. A1 in paper IV
-        
+
         x_spline = CubicSpline(s, x, bc_type='not-a-knot')
         dx_ds_spline = x_spline.derivative()
         d2x_ds2_spline = x_spline.derivative(2)
-       
+
         return s, x, x_spline, dx_ds_spline, d2x_ds2_spline
-    
+
     def integrate_forces(self, spline_points=1000, bc_type='not-a-knot'):
         """
         Use spline fit to integrate forces along MEP to approximate
@@ -854,7 +861,7 @@ class NEBOptimizer(Optimizer):
                  C1=1e-2,
                  C2=2.0):
 
-        Optimizer.__init__(self, atoms=neb, restart=restart, 
+        Optimizer.__init__(self, atoms=neb, restart=restart,
                            logfile=logfile, trajectory=trajectory,
                            master=master,
                            append_trajectory=append_trajectory,
@@ -937,26 +944,33 @@ class NEBOptimizer(Optimizer):
             except OptimizerConvergenceError:
                 return False
         elif self.method == 'krylov':
-            res = root(self.force_function,
-                       self.get_dofs(),
-                       method='krylov',
-                       options={'disp': True, 'fatol': fmax, 'maxiter': steps},
-                       callback=self.callback)
-            if res.success:
-                self.set_dofs(res.x)
-                return True
-            else:
-                return False
+            return self.thing2()
+        else:
+            return self.thing3()
+
+    def thing2(self):
+        res = root(self.force_function,
+                   self.get_dofs(),
+                   method='krylov',
+                   options={'disp': True, 'fatol': fmax, 'maxiter': steps},
+                   callback=self.callback)
+        if res.success:
+            self.set_dofs(res.x)
+            return True
         else:
-            X = self.get_dofs()
-            for step in range(steps):
-                F = self.neb.force_function(X)
-                if self.neb.get_residual() <= fmax:
-                    return True
-                X += self.alpha * F
-                self.callback(X)
             return False
-            
+
+
+    def thing3(self):
+        X = self.get_dofs()
+        for step in range(steps):
+            F = self.neb.force_function(X)
+            if self.neb.get_residual() <= fmax:
+                return True
+            X += self.alpha * F
+            self.callback(X)
+        return False
+
 
 
 class IDPP(Calculator):
diff --git a/ase/optimize/precon/precon.py b/ase/optimize/precon/precon.py
index 74d9927a5..b7b40e7b0 100644
--- a/ase/optimize/precon/precon.py
+++ b/ase/optimize/precon/precon.py
@@ -22,7 +22,7 @@ from ase.optimize.precon.neighbors import (get_neighbours,
 try:
     from pyamg import smoothed_aggregation_solver
     have_pyamg = True
-    
+
     def create_pyamg_solver(P, max_levels=15):
         return smoothed_aggregation_solver(
             P, B=None,
@@ -40,7 +40,7 @@ try:
             max_levels=max_levels,
             max_coarse=300,
             coarse_solver='pinv')
-    
+
 except ImportError:
     have_pyamg = False
 
@@ -60,13 +60,13 @@ class Precon(ABC):
 
         Args:
             atoms: the Atoms object used to create the preconditioner.
-            
+
             reinitialize: if True, parameters of the preconditioner
                 will be recalculated before the preconditioner matrix is
                 created. If False, they will be calculated only when they
                 do not currently have a value (ie, the first time this
                 function is called).
-            
+
 
         Returns:
             P: A sparse scipy csr_matrix. BE AWARE that using
@@ -82,7 +82,7 @@ class Precon(ABC):
         Return the result of applying P to a vector x
         """
         ...
-        
+
     def dot(self, x, y):
         """
         Return the preconditioned dot product <P x, y>
@@ -90,13 +90,13 @@ class Precon(ABC):
         Uses 128-bit floating point math for vector dot products
         """
         return longsum(self.Pdot(x) * y)
-    
+
     def norm(self, x):
         """
         Return the P-norm of x, where |x|_P = sqrt(<Px, x>)
         """
         return np.sqrt(self.dot(x, x))
-                
+
     def vdot(self, x, y):
         return self.dot(x.reshape(-1),
                         y.reshape(-1))
@@ -129,11 +129,11 @@ class Precon(ABC):
         residual = np.linalg.norm(forces, np.inf)
         precon_forces = self.solve(forces)
         return precon_forces, residual
-    
+
     @abstractmethod
     def copy(self):
         ...
-        
+
     @abstractmethod
     def asarray(self):
         """
@@ -259,11 +259,11 @@ class SparsePrecon(Precon):
             raise ValueError('Dimension must be at least 1')
         self.dim = dim
         self.logfile = Logfile(logfile)
-        
+
         if rng is None:
             rng = np.random.RandomState()
         self.rng = rng
-            
+
     def copy(self):
         return copy.deepcopy(self)
 
@@ -445,10 +445,10 @@ class SparsePrecon(Precon):
 
         self.P = None  # force a rebuild with new mu (there may be fixed atoms)
         return (self.mu, self.mu_c)
-    
+
     def asarray(self):
         return np.array(self.P.todense())
-    
+
     def one_dim_to_ndim(self, csc_P, N):
         """Expand an N x N precon matrix to self.dim*N x self.dim * N
 
@@ -576,7 +576,7 @@ class SparseCoeffPrecon(SparsePrecon):
 
         self.P = self.one_dim_to_ndim(csc_P, N)
         self.create_solver()
-    
+
     def make_precon(self, atoms, reinitialize=None):
         if self.r_NN is None:
             self.r_NN = estimate_nearest_neighbour_distance(atoms)
@@ -632,13 +632,7 @@ class SparseCoeffPrecon(SparsePrecon):
         ...
 
 
-class Pfrommer(Precon):
-    """Use initial guess for inverse Hessian from Pfrommer et al. as a
-    simple preconditioner
-
-    J. Comput. Phys. vol 131 p233-240 (1997)
-    """
-
+class PfrommerParameters:
     def __init__(self, bulk_modulus=500 * units.GPa, phonon_frequency=50 * THz,
                  apply_positions=True, apply_cell=True):
         """
@@ -648,12 +642,18 @@ class Pfrommer(Precon):
         self.phonon_frequency = phonon_frequency
         self.apply_positions = apply_positions
         self.apply_cell = apply_cell
-        self.H0 = None
+
+    def copy(self):
+        return PfrommerParameters(
+            self.bulk_modulus,
+            self.phonon_frequency,
+            self.apply_positions,
+            self.apply_cell)
 
     def make_precon(self, atoms, reinitialize=None):
-        if self.H0 is not None:
+        #if self.H0 is not None:
             # only build H0 on first call
-            return
+        #    return
 
         variable_cell = False
         if isinstance(atoms, Filter):
@@ -673,8 +673,27 @@ class Pfrommer(Precon):
                 coeff = 1.0 / (3 * self.bulk_modulus)
             blocks.append(np.diag([coeff] * 9))
 
-        self.H0 = sparse.block_diag(blocks, format='csr')
-        return
+        H0 = sparse.block_diag(blocks, format='csr')
+        return Pfrommer(H0, self)
+
+
+class Pfrommer(Precon):
+    """Use initial guess for inverse Hessian from Pfrommer et al. as a
+    simple preconditioner
+
+    J. Comput. Phys. vol 131 p233-240 (1997)
+    """
+
+    def __init__(self, bulk_modulus=500 * units.GPa, phonon_frequency=50 * THz,
+                 apply_positions=True, apply_cell=True):
+        """
+        Default bulk modulus is 500 GPa and default phonon frequency is 50 THz
+        """
+        #self.bulk_modulus = bulk_modulus
+        #self.phonon_frequency = phonon_frequency
+        #self.apply_positions = apply_positions
+        #self.apply_cell = apply_cell
+        self.H0 = H0
 
     def Pdot(self, x):
         return self.H0.solve(x)
@@ -682,13 +701,7 @@ class Pfrommer(Precon):
     def solve(self, x):
         y = self.H0.dot(x)
         return y
-    
-    def copy(self):
-        return Pfrommer(self.bulk_modulus,
-                        self.phonon_frequency,
-                        self.apply_positions,
-                        self.apply_cell)
-    
+
     def asarray(self):
         return np.array(self.H0.todense())
 
@@ -709,7 +722,7 @@ class IdentityPrecon(Precon):
 
     def copy(self):
         return IdentityPrecon()
-    
+
     def asarray(self):
         return np.eye(3 * len(self.atoms))
 
@@ -784,8 +797,8 @@ def ijk_to_x(i, j, k):
          3 * j, 3 * j + 1, 3 * j + 2,
          3 * k, 3 * k + 1, 3 * k + 2]
     return x
-    
-    
+
+
 def ijkl_to_x(i, j, k, l):
     x = [3 * i, 3 * i + 1, 3 * i + 2,
          3 * j, 3 * j + 1, 3 * j + 2,
@@ -839,7 +852,7 @@ class FF(SparsePrecon):
         self._make_sparse_precon(atoms, force_stab=self.force_stab)
         self.logfile.write('--- Precon created in %s seconds ---\n'
                            % (time.time() - start_time))
-        
+
     def add_morse(self, morse, atoms, row, col, data, conn=None):
         if self.hessian == 'reduced':
             i, j, Hx = ff.get_morse_potential_reduced_hessian(
@@ -856,7 +869,7 @@ class FF(SparsePrecon):
         if conn is not None:
             conn[i, j] = True
             conn[j, i] = True
-                        
+
     def add_bond(self, bond, atoms, row, col, data, conn=None):
         if self.hessian == 'reduced':
             i, j, Hx = ff.get_bond_potential_reduced_hessian(
@@ -873,7 +886,7 @@ class FF(SparsePrecon):
         if conn is not None:
             conn[i, j] = True
             conn[j, i] = True
-        
+
     def add_angle(self, angle, atoms, row, col, data, conn=None):
         if self.hessian == 'reduced':
             i, j, k, Hx = ff.get_angle_potential_reduced_hessian(
@@ -910,7 +923,7 @@ class FF(SparsePrecon):
                 j, k] = conn[j, l] = conn[k, l] = True
             conn[j, i] = conn[k, i] = conn[l, i] = conn[
                 k, j] = conn[l, j] = conn[l, k] = True
-        
+
     def _make_sparse_precon(self, atoms, initial_assembly=False,
                             force_stab=False):
         N = len(atoms)
@@ -922,7 +935,7 @@ class FF(SparsePrecon):
         if self.morses is not None:
             for morse in self.morses:
                 self.add_morse(morse, atoms, row, col, data)
-                
+
         if self.bonds is not None:
             for bond in self.bonds:
                 self.add_bond(bond, atoms, row, col, data)
@@ -1174,10 +1187,10 @@ def make_precon(precon, atoms=None, **kwargs):
     ----------
     precon - one of 'C1', 'Exp', 'Pfrommer', 'FF', 'Exp_FF', 'ID', None
              or an instance of a subclass of `ase.optimize.precon.Precon`
-             
+
     atoms - ase.atoms.Atoms instance, optional
             If present, build apreconditioner for this Atoms object
-            
+
     **kwargs - additional keyword arguments to pass to Precon constructor
 
     Returns
diff --git a/ase/utils/ff.py b/ase/utils/ff.py
index a5b823df0..5f4c5ac83 100644
--- a/ase/utils/ff.py
+++ b/ase/utils/ff.py
@@ -26,6 +26,20 @@ class Morse:
         self.r0 = r0
         self.r = None
 
+    def get_bond_potential_value(self, atoms):
+        i = self.atomi
+        j = self.atomj
+
+        rij = rel_pos_pbc(atoms, i, j)
+        dij = linalg.norm(rij)
+
+        v = 0.5*self.k*(dij-self.b0)**2
+
+        self.b = dij
+
+        return i, j, v
+
+
 class Bond:
     def __init__(self, atomi, atomj, k, b0, 
                  alpha=None, rref=None):
@@ -211,18 +225,6 @@ def get_morse_potential_reduced_hessian(atoms, morse):
 
     return i, j, Hx
 
-def get_bond_potential_value(atoms, bond):
-    i = bond.atomi
-    j = bond.atomj
-
-    rij = rel_pos_pbc(atoms, i, j)
-    dij = linalg.norm(rij)
-
-    v = 0.5*bond.k*(dij-bond.b0)**2
-
-    bond.b = dij
-
-    return i, j, v
 
 def get_bond_potential_gradient(atoms, bond):
     i = bond.atomi