Coverage for ase / optimize / precon / precon.py: 85.86%
665 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3"""
4Implementation of the Precon abstract base class and subclasses
5"""
7import copy
8import sys
9import time
10import warnings
11from abc import ABC, abstractmethod
13import numpy as np
14from scipy import sparse
15from scipy.interpolate import CubicSpline
16from scipy.sparse.linalg import spsolve
18import ase.units as units
19import ase.utils.ff as ff
20from ase.constraints import FixAtoms
21from ase.filters import Filter
22from ase.geometry import find_mic
23from ase.neighborlist import neighbor_list
24from ase.optimize.precon.neighbors import (
25 estimate_nearest_neighbour_distance,
26 get_neighbours,
27)
28from ase.utils import longsum, tokenize_version
30try:
31 import pyamg
32except ImportError:
33 have_pyamg = False
34else:
35 have_pyamg = True
38def create_pyamg_solver(P, max_levels=15):
39 if not have_pyamg:
40 raise RuntimeError(
41 'pyamg not available: install with `pip install pyamg`')
42 filter_key = 'filter'
43 if tokenize_version(pyamg.__version__) >= tokenize_version('4.2.1'):
44 filter_key = 'filter_entries'
45 return pyamg.smoothed_aggregation_solver(
46 P, B=None,
47 strength=('symmetric', {'theta': 0.0}),
48 smooth=(
49 'jacobi', {filter_key: True, 'weighting': 'local'}),
50 improve_candidates=([('block_gauss_seidel',
51 {'sweep': 'symmetric', 'iterations': 4})] +
52 [None] * (max_levels - 1)),
53 aggregate='standard',
54 presmoother=('block_gauss_seidel',
55 {'sweep': 'symmetric', 'iterations': 1}),
56 postsmoother=('block_gauss_seidel',
57 {'sweep': 'symmetric', 'iterations': 1}),
58 max_levels=max_levels,
59 max_coarse=300,
60 coarse_solver='pinv')
63THz = 1e12 * 1. / units.s
66class Precon(ABC):
68 @abstractmethod
69 def make_precon(self, atoms, reinitialize=None):
70 """
71 Create a preconditioner matrix based on the passed set of atoms.
73 Creates a general-purpose preconditioner for use with optimization
74 algorithms, based on examining distances between pairs of atoms in the
75 lattice. The matrix will be stored in the attribute self.P and
76 returned.
78 Args:
79 atoms: the Atoms object used to create the preconditioner.
81 reinitialize: if True, parameters of the preconditioner
82 will be recalculated before the preconditioner matrix is
83 created. If False, they will be calculated only when they
84 do not currently have a value (ie, the first time this
85 function is called).
87 Returns
88 -------
89 P: A sparse scipy csr_matrix. BE AWARE that using
90 numpy.dot() with sparse matrices will result in
91 errors/incorrect results - use the .dot method directly
92 on the matrix instead.
93 """
94 ...
96 @abstractmethod
97 def Pdot(self, x):
98 """
99 Return the result of applying P to a vector x
100 """
101 ...
103 def dot(self, x, y):
104 """
105 Return the preconditioned dot product <P x, y>
107 Uses 128-bit floating point math for vector dot products
108 """
109 return longsum(self.Pdot(x) * y)
111 def norm(self, x):
112 """
113 Return the P-norm of x, where |x|_P = sqrt(<Px, x>)
114 """
115 return np.sqrt(self.dot(x, x))
117 def vdot(self, x, y):
118 return self.dot(x.reshape(-1),
119 y.reshape(-1))
121 @abstractmethod
122 def solve(self, x):
123 """
124 Solve the (sparse) linear system P x = y and return y
125 """
126 ...
128 def apply(self, forces, atoms):
129 """
130 Convenience wrapper that combines make_precon() and solve()
132 Parameters
133 ----------
134 forces: array
135 (len(atoms)*3) array of input forces
136 atoms: ase.atoms.Atoms
138 Returns
139 -------
140 precon_forces: array
141 (len(atoms), 3) array of preconditioned forces
142 residual: float
143 inf-norm of original forces, i.e. maximum absolute force
144 """
145 self.make_precon(atoms)
146 residual = np.linalg.norm(forces, np.inf)
147 precon_forces = self.solve(forces)
148 return precon_forces, residual
150 @abstractmethod
151 def copy(self):
152 ...
154 @abstractmethod
155 def asarray(self):
156 """
157 Array representation of preconditioner, as a dense matrix
158 """
159 ...
162class Logfile:
163 def __init__(self, logfile=None):
164 if isinstance(logfile, str):
165 if logfile == "-":
166 logfile = sys.stdout
167 else:
168 logfile = open(logfile, "a")
169 self.logfile = logfile
171 def write(self, *args):
172 if self.logfile is None:
173 return
174 self.logfile.write(*args)
177class SparsePrecon(Precon):
178 def __init__(self, r_cut=None, r_NN=None,
179 mu=None, mu_c=None,
180 dim=3, c_stab=0.1, force_stab=False,
181 reinitialize=False, array_convention='C',
182 solver="auto", solve_tol=1e-8,
183 apply_positions=True, apply_cell=True,
184 estimate_mu_eigmode=False, logfile=None, rng=None,
185 neighbour_list=neighbor_list):
186 """Initialise a preconditioner object based on passed parameters.
188 Parameters
189 ----------
190 r_cut: float. This is a cut-off radius. The preconditioner matrix
191 will be created by considering pairs of atoms that are within a
192 distance r_cut of each other. For a regular lattice, this is
193 usually taken somewhere between the first- and second-nearest
194 neighbour distance. If r_cut is not provided, default is
195 2 * r_NN (see below)
196 r_NN: nearest neighbour distance. If not provided, this is
197 calculated
198 from input structure.
199 mu: float
200 energy scale for position degreees of freedom. If `None`, mu
201 is precomputed using finite difference derivatives.
202 mu_c: float
203 energy scale for cell degreees of freedom. Also precomputed
204 if None.
205 estimate_mu_eigmode:
206 If True, estimates mu based on the lowest eigenmodes of
207 unstabilised preconditioner. If False it uses the sine based
208 approach.
209 dim: int; dimensions of the problem
210 c_stab: float. The diagonal of the preconditioner matrix will have
211 a stabilisation constant added, which will be the value of
212 c_stab times mu.
213 force_stab:
214 If True, always add the stabilisation to diagnonal, regardless
215 of the presence of fixed atoms.
216 reinitialize: if True, the value of mu will be recalculated when
217 self.make_precon is called. This can be overridden in specific
218 cases with reinitialize argument in self.make_precon. If it
219 is set to True here, the value passed for mu will be
220 irrelevant unless reinitialize is set to False the first time
221 make_precon is called.
222 array_convention: Either 'C' or 'F' for Fortran; this will change
223 the preconditioner to reflect the ordering of the indices in
224 the vector it will operate on. The C convention assumes the
225 vector will be arranged atom-by-atom (ie [x1, y1, z1, x2, ...])
226 while the F convention assumes it will be arranged component
227 by component (ie [x1, x2, ..., y1, y2, ...]).
228 solver: One of "auto", "direct" or "pyamg", specifying whether to
229 use a direct sparse solver or PyAMG to solve P x = y.
230 Default is "auto" which uses PyAMG if available, falling
231 back to sparse solver if not. solve_tol: tolerance used for
232 PyAMG sparse linear solver, if available.
233 apply_positions: bool
234 if True, apply preconditioner to position DoF
235 apply_cell: bool
236 if True, apply preconditioner to cell DoF
237 logfile: file object or str
238 If *logfile* is a string, a file with that name will be opened.
239 Use '-' for stdout.
240 rng: None or np.random.RandomState instance
241 Random number generator to use for initialising pyamg solver
242 neighbor_list: function (optional). Optionally replace the built-in
243 ASE neighbour list with an alternative with the same call
244 signature, e.g. `matscipy.neighbours.neighbour_list`.
246 Raises
247 ------
248 ValueError for problem with arguments
250 """
251 self.r_NN = r_NN
252 self.r_cut = r_cut
253 self.mu = mu
254 self.mu_c = mu_c
255 self.estimate_mu_eigmode = estimate_mu_eigmode
256 self.c_stab = c_stab
257 self.force_stab = force_stab
258 self.array_convention = array_convention
259 self.reinitialize = reinitialize
260 self.P = None
261 self.old_positions = None
263 use_pyamg = False
264 if solver == "auto":
265 use_pyamg = have_pyamg
266 elif solver == "direct":
267 use_pyamg = False
268 elif solver == "pyamg":
269 if not have_pyamg:
270 raise RuntimeError("solver='pyamg', PyAMG can't be imported!")
271 use_pyamg = True
272 else:
273 raise ValueError('unknown solver - '
274 'should be "auto", "direct" or "pyamg"')
276 self.use_pyamg = use_pyamg
277 self.solve_tol = solve_tol
278 self.apply_positions = apply_positions
279 self.apply_cell = apply_cell
281 if dim < 1:
282 raise ValueError('Dimension must be at least 1')
283 self.dim = dim
284 self.logfile = Logfile(logfile)
286 if rng is None:
287 rng = np.random.RandomState()
288 self.rng = rng
290 self.neighbor_list = neighbor_list
292 def copy(self):
293 return copy.deepcopy(self)
295 def Pdot(self, x):
296 return self.P.dot(x)
298 def solve(self, x):
299 start_time = time.time()
300 if self.use_pyamg and have_pyamg:
301 y = self.ml.solve(x, x0=self.rng.random(self.P.shape[0]),
302 tol=self.solve_tol,
303 accel='cg',
304 maxiter=300,
305 cycle='W')
306 else:
307 y = spsolve(self.P, x)
308 self.logfile.write(
309 f'--- Precon applied in {(time.time() - start_time)} seconds ---\n')
310 return y
312 def estimate_mu(self, atoms, H=None):
313 r"""Estimate optimal preconditioner coefficient \mu
315 \mu is estimated from a numerical solution of
317 [dE(p+v) - dE(p)] \cdot v = \mu < P1 v, v >
319 with perturbation
321 v(x,y,z) = H P_lowest_nonzero_eigvec(x, y, z)
323 or
325 v(x,y,z) = H (sin(x / Lx), sin(y / Ly), sin(z / Lz))
327 After the optimal \mu is found, self.mu will be set to its value.
329 If `atoms` is an instance of Filter an additional \mu_c
330 will be computed for the cell degrees of freedom .
332 Args:
333 atoms: Atoms object for initial system
335 H: 3x3 array or None
336 Magnitude of deformation to apply.
337 Default is 1e-2*rNN*np.eye(3)
339 Returns
340 -------
341 mu : float
342 mu_c : float or None
343 """
344 logfile = self.logfile
346 if self.dim != 3:
347 raise ValueError('Automatic calculation of mu only possible for '
348 'three-dimensional preconditioners. Try setting '
349 'mu manually instead.')
351 if self.r_NN is None:
352 self.r_NN = estimate_nearest_neighbour_distance(atoms,
353 self.neighbor_list)
355 # deformation matrix, default is diagonal
356 if H is None:
357 H = 1e-2 * self.r_NN * np.eye(3)
359 # compute perturbation
360 p = atoms.get_positions()
362 if self.estimate_mu_eigmode:
363 self.mu = 1.0
364 self.mu_c = 1.0
365 c_stab = self.c_stab
366 self.c_stab = 0.0
368 if isinstance(atoms, Filter):
369 n = len(atoms.atoms)
370 else:
371 n = len(atoms)
372 self._make_sparse_precon(atoms, initial_assembly=True)
373 self.P = self.P[:3 * n, :3 * n]
374 eigvals, eigvecs = sparse.linalg.eigsh(self.P, k=4, which='SM')
376 logfile.write('estimate_mu(): lowest 4 eigvals = %f %f %f %f\n' %
377 (eigvals[0], eigvals[1], eigvals[2], eigvals[3]))
378 # check eigenvalues
379 if any(eigvals[0:3] > 1e-6):
380 raise ValueError('First 3 eigenvalues of preconditioner matrix'
381 'do not correspond to translational modes.')
382 elif eigvals[3] < 1e-6:
383 raise ValueError('Fourth smallest eigenvalue of '
384 'preconditioner matrix '
385 'is too small, increase r_cut.')
387 x = np.zeros(n)
388 for i in range(n):
389 x[i] = eigvecs[:, 3][3 * i]
390 x = x / np.linalg.norm(x)
391 if x[0] < 0:
392 x = -x
394 v = np.zeros(3 * len(atoms))
395 for i in range(n):
396 v[3 * i] = x[i]
397 v[3 * i + 1] = x[i]
398 v[3 * i + 2] = x[i]
399 v = v / np.linalg.norm(v)
400 v = v.reshape((-1, 3))
402 self.c_stab = c_stab
403 else:
404 Lx, Ly, Lz = (p[:, i].max() - p[:, i].min() for i in range(3))
405 logfile.write('estimate_mu(): Lx=%.1f Ly=%.1f Lz=%.1f\n' %
406 (Lx, Ly, Lz))
408 x, y, z = p.T
409 # sine_vr = [np.sin(x/Lx), np.sin(y/Ly), np.sin(z/Lz)], but we need
410 # to take into account the possibility that one of Lx/Ly/Lz is
411 # zero.
412 sine_vr = [x, y, z]
414 for i, L in enumerate([Lx, Ly, Lz]):
415 if L == 0:
416 warnings.warn(
417 'Cell length L[%d] == 0. Setting H[%d,%d] = 0.' %
418 (i, i, i))
419 H[i, i] = 0.0
420 else:
421 sine_vr[i] = np.sin(sine_vr[i] / L)
423 v = np.dot(H, sine_vr).T
425 natoms = len(atoms)
426 if isinstance(atoms, Filter):
427 natoms = len(atoms.atoms)
428 eps = H / self.r_NN
429 v[natoms:, :] = eps
431 v1 = v.reshape(-1)
433 # compute LHS
434 dE_p = -atoms.get_forces().reshape(-1)
435 atoms_v = atoms.copy()
436 atoms_v.calc = atoms.calc
437 if isinstance(atoms, Filter):
438 atoms_v = atoms.__class__(atoms_v)
439 if hasattr(atoms, 'constant_volume'):
440 atoms_v.constant_volume = atoms.constant_volume
441 atoms_v.set_positions(p + v)
442 dE_p_plus_v = -atoms_v.get_forces().reshape(-1)
444 # compute left hand side
445 LHS = (dE_p_plus_v - dE_p) * v1
447 # assemble P with \mu = 1
448 self.mu = 1.0
449 self.mu_c = 1.0
451 self._make_sparse_precon(atoms, initial_assembly=True)
453 # compute right hand side
454 RHS = self.P.dot(v1) * v1
456 # use partial sums to compute separate mu for positions and cell DoFs
457 self.mu = longsum(LHS[:3 * natoms]) / longsum(RHS[:3 * natoms])
458 if self.mu < 1.0:
459 logfile.write('estimate_mu(): mu (%.3f) < 1.0, '
460 'capping at mu=1.0' % self.mu)
461 self.mu = 1.0
463 if isinstance(atoms, Filter):
464 self.mu_c = longsum(LHS[3 * natoms:]) / longsum(RHS[3 * natoms:])
465 if self.mu_c < 1.0:
466 logfile.write('estimate_mu(): mu_c (%.3f) < 1.0, '
467 'capping at mu_c=1.0\n' % self.mu_c)
468 self.mu_c = 1.0
470 logfile.write('estimate_mu(): mu=%r, mu_c=%r\n' %
471 (self.mu, self.mu_c))
473 self.P = None # force a rebuild with new mu (there may be fixed atoms)
474 return (self.mu, self.mu_c)
476 def asarray(self):
477 return np.array(self.P.todense())
479 def one_dim_to_ndim(self, csc_P, N):
480 """
481 Expand an N x N precon matrix to self.dim*N x self.dim * N
483 Args:
484 csc_P (sparse matrix): N x N sparse matrix, in CSC format
485 """
486 start_time = time.time()
487 if self.dim == 1:
488 P = csc_P
489 elif self.array_convention == 'F':
490 csc_P = csc_P.tocsr()
491 P = csc_P
492 for _ in range(self.dim - 1):
493 P = sparse.block_diag((P, csc_P)).tocsr()
494 else:
495 # convert back to triplet and read the arrays
496 csc_P = csc_P.tocoo()
497 i = csc_P.row * self.dim
498 j = csc_P.col * self.dim
499 z = csc_P.data
501 # N-dimensionalise, interlaced coordinates
502 I = np.hstack([i + d for d in range(self.dim)])
503 J = np.hstack([j + d for d in range(self.dim)])
504 Z = np.hstack([z for _ in range(self.dim)])
505 P = sparse.csc_matrix((Z, (I, J)),
506 shape=(self.dim * N, self.dim * N))
507 P = P.tocsr()
508 self.logfile.write(
509 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
510 return P
512 def create_solver(self):
513 if self.use_pyamg and have_pyamg:
514 start_time = time.time()
515 self.ml = create_pyamg_solver(self.P)
516 self.logfile.write(
517 f'--- multi grid solver created in {(time.time() - start_time)}'
518 ' ---\n')
521class SparseCoeffPrecon(SparsePrecon):
522 def _make_sparse_precon(self, atoms, initial_assembly=False,
523 force_stab=False):
524 """Create a sparse preconditioner matrix based on the passed atoms.
526 Creates a general-purpose preconditioner for use with optimization
527 algorithms, based on examining distances between pairs of atoms in the
528 lattice. The matrix will be stored in the attribute self.P and
529 returned. Note that this function will use self.mu, whatever it is.
531 Args:
532 atoms: the Atoms object used to create the preconditioner.
534 Returns
535 -------
536 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
537 (where N is the number of atoms, and d is the value of self.dim).
538 BE AWARE that using numpy.dot() with this object will result in
539 errors/incorrect results - use the .dot method directly on the
540 sparse matrix instead.
542 """
543 logfile = self.logfile
544 logfile.write('creating sparse precon: initial_assembly=%r, '
545 'force_stab=%r, apply_positions=%r, apply_cell=%r\n' %
546 (initial_assembly, force_stab, self.apply_positions,
547 self.apply_cell))
549 N = len(atoms)
550 diag_i = np.arange(N, dtype=int)
551 start_time = time.time()
552 if self.apply_positions:
553 # compute neighbour list
554 i, j, rij, fixed_atoms = get_neighbours(
555 atoms, self.r_cut,
556 neighbor_list=self.neighbor_list)
557 logfile.write(
558 f'--- neighbour list created in {(time.time() - start_time)} s '
559 '--- \n')
561 # compute entries in triplet format: without the constraints
562 start_time = time.time()
563 coeff = self.get_coeff(rij)
564 diag_coeff = np.bincount(i, -coeff, minlength=N).astype(np.float64)
565 if force_stab or len(fixed_atoms) == 0:
566 logfile.write('adding stabilisation to precon')
567 diag_coeff += self.mu * self.c_stab
568 else:
569 diag_coeff = np.ones(N)
571 # precon is mu_c * identity for cell DoF
572 if isinstance(atoms, Filter):
573 if self.apply_cell:
574 diag_coeff[-3:] = self.mu_c
575 else:
576 diag_coeff[-3:] = 1.0
577 logfile.write(
578 f'--- computed triplet format in {(time.time() - start_time)} s '
579 '---\n')
581 if self.apply_positions and not initial_assembly:
582 # apply the constraints
583 start_time = time.time()
584 mask = np.ones(N)
585 mask[fixed_atoms] = 0.0
586 coeff *= mask[i] * mask[j]
587 diag_coeff[fixed_atoms] = 1.0
588 logfile.write(
589 f'--- applied fixed_atoms in {(time.time() - start_time)} s '
590 '---\n')
592 if self.apply_positions:
593 # remove zeros
594 start_time = time.time()
595 inz = np.nonzero(coeff)
596 i = np.hstack((i[inz], diag_i))
597 j = np.hstack((j[inz], diag_i))
598 coeff = np.hstack((coeff[inz], diag_coeff))
599 logfile.write(
600 f'--- remove zeros in {(time.time() - start_time)} s '
601 '---\n')
602 else:
603 i = diag_i
604 j = diag_i
605 coeff = diag_coeff
607 # create an N x N precon matrix in compressed sparse column (CSC) format
608 start_time = time.time()
609 csc_P = sparse.csc_matrix((coeff, (i, j)), shape=(N, N))
610 logfile.write(
611 f'--- created CSC matrix in {(time.time() - start_time)} s '
612 '---\n')
614 self.P = self.one_dim_to_ndim(csc_P, N)
615 self.create_solver()
617 def make_precon(self, atoms, reinitialize=None):
618 if self.r_NN is None:
619 self.r_NN = estimate_nearest_neighbour_distance(atoms,
620 self.neighbor_list)
622 if self.r_cut is None:
623 # This is the first time this function has been called, and no
624 # cutoff radius has been specified, so calculate it automatically.
625 self.r_cut = 2.0 * self.r_NN
626 elif self.r_cut < self.r_NN:
627 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
628 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
629 self.r_NN,
630 1.1 * self.r_NN))
631 warnings.warn(warning)
632 self.r_cut = 1.1 * self.r_NN
634 if reinitialize is None:
635 # The caller has not specified whether or not to recalculate mu,
636 # so the Precon's setting is used.
637 reinitialize = self.reinitialize
639 if self.mu is None:
640 # Regardless of what the caller has specified, if we don't
641 # currently have a value of mu, then we need one.
642 reinitialize = True
644 if reinitialize:
645 self.estimate_mu(atoms)
647 if self.P is not None:
648 real_atoms = atoms
649 if isinstance(atoms, Filter):
650 real_atoms = atoms.atoms
651 if self.old_positions is None:
652 self.old_positions = real_atoms.positions
653 displacement, _ = find_mic(real_atoms.positions -
654 self.old_positions,
655 real_atoms.cell, real_atoms.pbc)
656 self.old_positions = real_atoms.get_positions()
657 max_abs_displacement = abs(displacement).max()
658 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
659 (max_abs_displacement,
660 max_abs_displacement / self.r_NN))
661 if max_abs_displacement < 0.5 * self.r_NN:
662 return
664 start_time = time.time()
665 self._make_sparse_precon(atoms, force_stab=self.force_stab)
666 self.logfile.write(
667 f'--- Precon created in {(time.time() - start_time)} seconds '
668 '--- \n')
670 @abstractmethod
671 def get_coeff(self, r):
672 ...
675class Pfrommer(Precon):
676 """
677 Use initial guess for inverse Hessian from Pfrommer et al. as a
678 simple preconditioner
680 J. Comput. Phys. vol 131 p233-240 (1997)
681 """
683 def __init__(self, bulk_modulus=500 * units.GPa, phonon_frequency=50 * THz,
684 apply_positions=True, apply_cell=True):
685 """
686 Default bulk modulus is 500 GPa and default phonon frequency is 50 THz
687 """
688 self.bulk_modulus = bulk_modulus
689 self.phonon_frequency = phonon_frequency
690 self.apply_positions = apply_positions
691 self.apply_cell = apply_cell
692 self.H0 = None
694 def make_precon(self, atoms, reinitialize=None):
695 if self.H0 is not None:
696 # only build H0 on first call
697 return
699 variable_cell = False
700 if isinstance(atoms, Filter):
701 variable_cell = True
702 atoms = atoms.atoms
704 # position DoF
705 omega = self.phonon_frequency
706 mass = atoms.get_masses().mean()
707 block = np.eye(3) / (mass * omega**2)
708 blocks = [block] * len(atoms)
710 # cell DoF
711 if variable_cell:
712 coeff = 1.0
713 if self.apply_cell:
714 coeff = 1.0 / (3 * self.bulk_modulus)
715 blocks.append(np.diag([coeff] * 9))
717 self.H0 = sparse.block_diag(blocks, format='csr')
718 return
720 def Pdot(self, x):
721 return self.H0.solve(x)
723 def solve(self, x):
724 y = self.H0.dot(x)
725 return y
727 def copy(self):
728 return Pfrommer(self.bulk_modulus,
729 self.phonon_frequency,
730 self.apply_positions,
731 self.apply_cell)
733 def asarray(self):
734 return np.array(self.H0.todense())
737class IdentityPrecon(Precon):
738 """
739 Dummy preconditioner which does not modify forces
740 """
742 def make_precon(self, atoms, reinitialize=None):
743 self.atoms = atoms
745 def Pdot(self, x):
746 return x
748 def solve(self, x):
749 return x
751 def copy(self):
752 return IdentityPrecon()
754 def asarray(self):
755 return np.eye(3 * len(self.atoms))
758class C1(SparseCoeffPrecon):
759 """Creates matrix by inserting a constant whenever r_ij is less than r_cut.
760 """
762 def __init__(self, r_cut=None, mu=None, mu_c=None, dim=3, c_stab=0.1,
763 force_stab=False,
764 reinitialize=False, array_convention='C',
765 solver="auto", solve_tol=1e-9,
766 apply_positions=True, apply_cell=True, logfile=None):
767 super().__init__(r_cut=r_cut, mu=mu, mu_c=mu_c,
768 dim=dim, c_stab=c_stab,
769 force_stab=force_stab,
770 reinitialize=reinitialize,
771 array_convention=array_convention,
772 solver=solver, solve_tol=solve_tol,
773 apply_positions=apply_positions,
774 apply_cell=apply_cell,
775 logfile=logfile)
777 def get_coeff(self, r):
778 return -self.mu * np.ones_like(r)
781class Exp(SparseCoeffPrecon):
782 """Creates matrix with values decreasing exponentially with distance.
783 """
785 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
786 dim=3, c_stab=0.1,
787 force_stab=False, reinitialize=False, array_convention='C',
788 solver="auto", solve_tol=1e-9,
789 apply_positions=True, apply_cell=True,
790 estimate_mu_eigmode=False, logfile=None):
791 """
792 Initialise an Exp preconditioner with given parameters.
794 Args:
795 r_cut, mu, c_stab, dim, sparse, reinitialize, array_convention: see
796 precon.__init__()
797 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
798 """
799 super().__init__(r_cut=r_cut, r_NN=r_NN,
800 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
801 force_stab=force_stab,
802 reinitialize=reinitialize,
803 array_convention=array_convention,
804 solver=solver,
805 solve_tol=solve_tol,
806 apply_positions=apply_positions,
807 apply_cell=apply_cell,
808 estimate_mu_eigmode=estimate_mu_eigmode,
809 logfile=logfile)
811 self.A = A
813 def get_coeff(self, r):
814 return -self.mu * np.exp(-self.A * (r / self.r_NN - 1))
817def ij_to_x(i, j):
818 x = [3 * i, 3 * i + 1, 3 * i + 2,
819 3 * j, 3 * j + 1, 3 * j + 2]
820 return x
823def ijk_to_x(i, j, k):
824 x = [3 * i, 3 * i + 1, 3 * i + 2,
825 3 * j, 3 * j + 1, 3 * j + 2,
826 3 * k, 3 * k + 1, 3 * k + 2]
827 return x
830def ijkl_to_x(i, j, k, l):
831 x = [3 * i, 3 * i + 1, 3 * i + 2,
832 3 * j, 3 * j + 1, 3 * j + 2,
833 3 * k, 3 * k + 1, 3 * k + 2,
834 3 * l, 3 * l + 1, 3 * l + 2]
835 return x
838def apply_fixed(atoms, P):
839 fixed_atoms = []
840 for constraint in atoms.constraints:
841 if isinstance(constraint, FixAtoms):
842 fixed_atoms.extend(list(constraint.index))
843 else:
844 raise TypeError(
845 'only FixAtoms constraints are supported by Precon class')
846 if len(fixed_atoms) != 0:
847 P = P.tolil()
848 for i in fixed_atoms:
849 P[i, :] = 0.0
850 P[:, i] = 0.0
851 P[i, i] = 1.0
852 return P
855class FF(SparsePrecon):
856 """Creates matrix using morse/bond/angle/dihedral force field parameters.
857 """
859 def __init__(self, dim=3, c_stab=0.1, force_stab=False,
860 array_convention='C', solver="auto", solve_tol=1e-9,
861 apply_positions=True, apply_cell=True,
862 hessian='spectral', morses=None, bonds=None, angles=None,
863 dihedrals=None, logfile=None):
864 """Initialise an FF preconditioner with given parameters.
866 Args:
867 dim, c_stab, force_stab, array_convention, use_pyamg, solve_tol:
868 see SparsePrecon.__init__()
869 morses: Morse instance
870 bonds: Bond instance
871 angles: Angle instance
872 dihedrals: Dihedral instance
873 """
875 if (morses is None and bonds is None and angles is None and
876 dihedrals is None):
877 raise ImportError(
878 'At least one of morses, bonds, angles or dihedrals must be '
879 'defined!')
881 super().__init__(dim=dim, c_stab=c_stab,
882 force_stab=force_stab,
883 array_convention=array_convention,
884 solver=solver,
885 solve_tol=solve_tol,
886 apply_positions=apply_positions,
887 apply_cell=apply_cell, logfile=logfile)
889 self.hessian = hessian
890 self.morses = morses
891 self.bonds = bonds
892 self.angles = angles
893 self.dihedrals = dihedrals
895 def make_precon(self, atoms, reinitialize=None):
896 start_time = time.time()
897 self._make_sparse_precon(atoms, force_stab=self.force_stab)
898 self.logfile.write(
899 f'--- Precon created in {(time.time() - start_time)} seconds '
900 '---\n')
902 def add_morse(self, morse, atoms, row, col, data, conn=None):
903 if self.hessian == 'reduced':
904 i, j, Hx = ff.get_morse_potential_reduced_hessian(
905 atoms, morse)
906 elif self.hessian == 'spectral':
907 i, j, Hx = ff.get_morse_potential_hessian(
908 atoms, morse, spectral=True)
909 else:
910 raise NotImplementedError('Not implemented hessian')
911 x = ij_to_x(i, j)
912 row.extend(np.repeat(x, 6))
913 col.extend(np.tile(x, 6))
914 data.extend(Hx.flatten())
915 if conn is not None:
916 conn[i, j] = True
917 conn[j, i] = True
919 def add_bond(self, bond, atoms, row, col, data, conn=None):
920 if self.hessian == 'reduced':
921 i, j, Hx = ff.get_bond_potential_reduced_hessian(
922 atoms, bond, self.morses)
923 elif self.hessian == 'spectral':
924 i, j, Hx = ff.get_bond_potential_hessian(
925 atoms, bond, self.morses, spectral=True)
926 else:
927 raise NotImplementedError('Not implemented hessian')
928 x = ij_to_x(i, j)
929 row.extend(np.repeat(x, 6))
930 col.extend(np.tile(x, 6))
931 data.extend(Hx.flatten())
932 if conn is not None:
933 conn[i, j] = True
934 conn[j, i] = True
936 def add_angle(self, angle, atoms, row, col, data, conn=None):
937 if self.hessian == 'reduced':
938 i, j, k, Hx = ff.get_angle_potential_reduced_hessian(
939 atoms, angle, self.morses)
940 elif self.hessian == 'spectral':
941 i, j, k, Hx = ff.get_angle_potential_hessian(
942 atoms, angle, self.morses, spectral=True)
943 else:
944 raise NotImplementedError('Not implemented hessian')
945 x = ijk_to_x(i, j, k)
946 row.extend(np.repeat(x, 9))
947 col.extend(np.tile(x, 9))
948 data.extend(Hx.flatten())
949 if conn is not None:
950 conn[i, j] = conn[i, k] = conn[j, k] = True
951 conn[j, i] = conn[k, i] = conn[k, j] = True
953 def add_dihedral(self, dihedral, atoms, row, col, data, conn=None):
954 if self.hessian == 'reduced':
955 i, j, k, l, Hx = \
956 ff.get_dihedral_potential_reduced_hessian(
957 atoms, dihedral, self.morses)
958 elif self.hessian == 'spectral':
959 i, j, k, l, Hx = ff.get_dihedral_potential_hessian(
960 atoms, dihedral, self.morses, spectral=True)
961 else:
962 raise NotImplementedError('Not implemented hessian')
963 x = ijkl_to_x(i, j, k, l)
964 row.extend(np.repeat(x, 12))
965 col.extend(np.tile(x, 12))
966 data.extend(Hx.flatten())
967 if conn is not None:
968 conn[i, j] = conn[i, k] = conn[i, l] = conn[
969 j, k] = conn[j, l] = conn[k, l] = True
970 conn[j, i] = conn[k, i] = conn[l, i] = conn[
971 k, j] = conn[l, j] = conn[l, k] = True
973 def _make_sparse_precon(self, atoms, initial_assembly=False,
974 force_stab=False):
975 N = len(atoms)
977 row = []
978 col = []
979 data = []
981 if self.morses is not None:
982 for morse in self.morses:
983 self.add_morse(morse, atoms, row, col, data)
985 if self.bonds is not None:
986 for bond in self.bonds:
987 self.add_bond(bond, atoms, row, col, data)
989 if self.angles is not None:
990 for angle in self.angles:
991 self.add_angle(angle, atoms, row, col, data)
993 if self.dihedrals is not None:
994 for dihedral in self.dihedrals:
995 self.add_dihedral(dihedral, atoms, row, col, data)
997 # add the diagonal
998 row.extend(range(self.dim * N))
999 col.extend(range(self.dim * N))
1000 data.extend([self.c_stab] * self.dim * N)
1002 # create the matrix
1003 start_time = time.time()
1004 self.P = sparse.csc_matrix(
1005 (data, (row, col)), shape=(self.dim * N, self.dim * N))
1006 self.logfile.write(
1007 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1009 self.P = apply_fixed(atoms, self.P)
1010 self.P = self.P.tocsr()
1011 self.logfile.write(
1012 f'--- N-dim precon created in {(time.time() - start_time)} s ---\n')
1013 self.create_solver()
1016class Exp_FF(Exp, FF):
1017 """Creates matrix with values decreasing exponentially with distance.
1018 """
1020 def __init__(self, A=3.0, r_cut=None, r_NN=None, mu=None, mu_c=None,
1021 dim=3, c_stab=0.1,
1022 force_stab=False, reinitialize=False, array_convention='C',
1023 solver="auto", solve_tol=1e-9,
1024 apply_positions=True, apply_cell=True,
1025 estimate_mu_eigmode=False,
1026 hessian='spectral', morses=None, bonds=None, angles=None,
1027 dihedrals=None, logfile=None):
1028 """Initialise an Exp+FF preconditioner with given parameters.
1030 Args:
1031 r_cut, mu, c_stab, dim, reinitialize, array_convention: see
1032 precon.__init__()
1033 A: coefficient in exp(-A*r/r_NN). Default is A=3.0.
1034 """
1035 if (morses is None and bonds is None and angles is None and
1036 dihedrals is None):
1037 raise ImportError(
1038 'At least one of morses, bonds, angles or dihedrals must '
1039 'be defined!')
1041 SparsePrecon.__init__(self, r_cut=r_cut, r_NN=r_NN,
1042 mu=mu, mu_c=mu_c, dim=dim, c_stab=c_stab,
1043 force_stab=force_stab,
1044 reinitialize=reinitialize,
1045 array_convention=array_convention,
1046 solver=solver,
1047 solve_tol=solve_tol,
1048 apply_positions=apply_positions,
1049 apply_cell=apply_cell,
1050 estimate_mu_eigmode=estimate_mu_eigmode,
1051 logfile=logfile)
1053 self.A = A
1054 self.hessian = hessian
1055 self.morses = morses
1056 self.bonds = bonds
1057 self.angles = angles
1058 self.dihedrals = dihedrals
1060 def make_precon(self, atoms, reinitialize=None):
1061 if self.r_NN is None:
1062 self.r_NN = estimate_nearest_neighbour_distance(atoms,
1063 self.neighbor_list)
1065 if self.r_cut is None:
1066 # This is the first time this function has been called, and no
1067 # cutoff radius has been specified, so calculate it automatically.
1068 self.r_cut = 2.0 * self.r_NN
1069 elif self.r_cut < self.r_NN:
1070 warning = ('WARNING: r_cut (%.2f) < r_NN (%.2f), '
1071 'increasing to 1.1*r_NN = %.2f' % (self.r_cut,
1072 self.r_NN,
1073 1.1 * self.r_NN))
1074 warnings.warn(warning)
1075 self.r_cut = 1.1 * self.r_NN
1077 if reinitialize is None:
1078 # The caller has not specified whether or not to recalculate mu,
1079 # so the Precon's setting is used.
1080 reinitialize = self.reinitialize
1082 if self.mu is None:
1083 # Regardless of what the caller has specified, if we don't
1084 # currently have a value of mu, then we need one.
1085 reinitialize = True
1087 if reinitialize:
1088 self.estimate_mu(atoms)
1090 if self.P is not None:
1091 real_atoms = atoms
1092 if isinstance(atoms, Filter):
1093 real_atoms = atoms.atoms
1094 if self.old_positions is None:
1095 self.old_positions = real_atoms.positions
1096 displacement, _ = find_mic(real_atoms.positions -
1097 self.old_positions,
1098 real_atoms.cell, real_atoms.pbc)
1099 self.old_positions = real_atoms.get_positions()
1100 max_abs_displacement = abs(displacement).max()
1101 self.logfile.write('max(abs(displacements)) = %.2f A (%.2f r_NN)' %
1102 (max_abs_displacement,
1103 max_abs_displacement / self.r_NN))
1104 if max_abs_displacement < 0.5 * self.r_NN:
1105 return
1107 # Create the preconditioner:
1108 start_time = time.time()
1109 self._make_sparse_precon(atoms, force_stab=self.force_stab)
1110 self.logfile.write(
1111 f'--- Precon created in {(time.time() - start_time)} seconds ---\n')
1113 def _make_sparse_precon(self, atoms, initial_assembly=False,
1114 force_stab=False):
1115 """Create a sparse preconditioner matrix based on the passed atoms.
1117 Args:
1118 atoms: the Atoms object used to create the preconditioner.
1120 Returns
1121 -------
1122 A scipy.sparse.csr_matrix object, representing a d*N by d*N matrix
1123 (where N is the number of atoms, and d is the value of self.dim).
1124 BE AWARE that using numpy.dot() with this object will result in
1125 errors/incorrect results - use the .dot method directly on the
1126 sparse matrix instead.
1128 """
1129 self.logfile.write('creating sparse precon: initial_assembly=%r, '
1130 'force_stab=%r, apply_positions=%r, '
1131 'apply_cell=%r\n' %
1132 (initial_assembly, force_stab,
1133 self.apply_positions, self.apply_cell))
1135 N = len(atoms)
1136 start_time = time.time()
1137 if self.apply_positions:
1138 # compute neighbour list
1139 i_list, j_list, rij_list, _fixed_atoms = get_neighbours(
1140 atoms, self.r_cut, self.neighbor_list)
1141 self.logfile.write(
1142 f'--- neighbour list created in {(time.time() - start_time)} s '
1143 '---\n')
1145 row = []
1146 col = []
1147 data = []
1149 # precon is mu_c*identity for cell DoF
1150 start_time = time.time()
1151 if isinstance(atoms, Filter):
1152 i = N - 3
1153 j = N - 2
1154 k = N - 1
1155 x = ijk_to_x(i, j, k)
1156 row.extend(x)
1157 col.extend(x)
1158 if self.apply_cell:
1159 data.extend(np.repeat(self.mu_c, 9))
1160 else:
1161 data.extend(np.repeat(self.mu_c, 9))
1162 self.logfile.write(
1163 f'--- computed triplet format in {(time.time() - start_time)} s '
1164 '---\n')
1166 conn = sparse.lil_matrix((N, N), dtype=bool)
1168 if self.apply_positions and not initial_assembly:
1169 if self.morses is not None:
1170 for morse in self.morses:
1171 self.add_morse(morse, atoms, row, col, data, conn)
1173 if self.bonds is not None:
1174 for bond in self.bonds:
1175 self.add_bond(bond, atoms, row, col, data, conn)
1177 if self.angles is not None:
1178 for angle in self.angles:
1179 self.add_angle(angle, atoms, row, col, data, conn)
1181 if self.dihedrals is not None:
1182 for dihedral in self.dihedrals:
1183 self.add_dihedral(dihedral, atoms, row, col, data, conn)
1185 if self.apply_positions:
1186 for i, j, rij in zip(i_list, j_list, rij_list):
1187 if not conn[i, j]:
1188 coeff = self.get_coeff(rij)
1189 x = ij_to_x(i, j)
1190 row.extend(x)
1191 col.extend(x)
1192 data.extend(3 * [-coeff] + 3 * [coeff])
1194 row.extend(range(self.dim * N))
1195 col.extend(range(self.dim * N))
1196 if initial_assembly:
1197 data.extend([self.mu * self.c_stab] * self.dim * N)
1198 else:
1199 data.extend([self.c_stab] * self.dim * N)
1201 # create the matrix
1202 start_time = time.time()
1203 self.P = sparse.csc_matrix(
1204 (data, (row, col)), shape=(self.dim * N, self.dim * N))
1205 self.logfile.write(
1206 f'--- created CSC matrix in {(time.time() - start_time)} s ---\n')
1208 if not initial_assembly:
1209 self.P = apply_fixed(atoms, self.P)
1211 self.P = self.P.tocsr()
1212 self.create_solver()
1215def make_precon(precon, atoms=None, **kwargs):
1216 """
1217 Construct preconditioner from a string and optionally build for Atoms
1219 Parameters
1220 ----------
1221 precon - one of 'C1', 'Exp', 'Pfrommer', 'FF', 'Exp_FF', 'ID', None
1222 or an instance of a subclass of `ase.optimize.precon.Precon`
1224 atoms - ase.atoms.Atoms instance, optional
1225 If present, build apreconditioner for this Atoms object
1227 **kwargs - additional keyword arguments to pass to Precon constructor
1229 Returns
1230 -------
1231 precon - instance of relevant subclass of `ase.optimize.precon.Precon`
1232 """
1233 lookup = {
1234 'C1': C1,
1235 'Exp': Exp,
1236 'Pfrommer': Pfrommer,
1237 'FF': FF,
1238 'Exp_FF': Exp_FF,
1239 'ID': IdentityPrecon,
1240 'IdentityPrecon': IdentityPrecon,
1241 None: IdentityPrecon
1242 }
1243 if isinstance(precon, str) or precon is None:
1244 cls = lookup[precon]
1245 precon = cls(**kwargs)
1246 if atoms is not None:
1247 precon.make_precon(atoms)
1248 return precon
1251class SplineFit:
1252 """
1253 Fit a cubic spline fit to images
1254 """
1256 def __init__(self, s, x):
1257 self._s = s
1258 self._x_data = x
1259 self._x = CubicSpline(self._s, x, bc_type='not-a-knot')
1260 self._dx_ds = self._x.derivative()
1261 self._d2x_ds2 = self._x.derivative(2)
1263 @property
1264 def s(self):
1265 return self._s
1267 @property
1268 def x_data(self):
1269 return self._x_data
1271 @property
1272 def x(self):
1273 return self._x
1275 @property
1276 def dx_ds(self):
1277 return self._dx_ds
1279 @property
1280 def d2x_ds2(self):
1281 return self._d2x_ds2
1284class PreconImages:
1285 def __init__(self, precon, images, **kwargs):
1286 """
1287 Wrapper for a list of Precon objects and associated images
1289 This is used when preconditioning a NEB object. Equation references
1290 refer to Paper IV in the :class:`ase.mep.NEB` documentation, i.e.
1292 S. Makri, C. Ortner and J. R. Kermode, J. Chem. Phys.
1293 150, 094109 (2019)
1294 https://dx.doi.org/10.1063/1.5064465
1296 Args:
1297 precon (str or list): preconditioner(s) to use
1298 images (list of Atoms): Atoms objects that define the state
1300 """
1301 self.images = images
1302 self._spline = None
1304 if isinstance(precon, list):
1305 if len(precon) != len(images):
1306 raise ValueError(f'length mismatch: len(precon)={len(precon)} '
1307 f'!= len(images)={len(images)}')
1308 self.precon = precon
1309 return
1310 P0 = make_precon(precon, images[0], **kwargs)
1311 self.precon = [P0]
1312 for image in images[1:]:
1313 P = P0.copy()
1314 P.make_precon(image)
1315 self.precon.append(P)
1317 def __len__(self):
1318 return len(self.precon)
1320 def __iter__(self):
1321 return iter(self.precon)
1323 def __getitem__(self, index):
1324 return self.precon[index]
1326 def apply(self, all_forces, index=None):
1327 """Apply preconditioners to stored images
1329 Args:
1330 all_forces (array): forces on images, shape (nimages, natoms, 3)
1331 index (slice, optional): Which images to include. Defaults to all.
1333 Returns
1334 -------
1335 precon_forces: array of preconditioned forces
1336 """
1337 if index is None:
1338 index = slice(None)
1339 precon_forces = []
1340 for precon, image, forces in zip(self.precon[index],
1341 self.images[index],
1342 all_forces):
1343 f_vec = forces.reshape(-1)
1344 pf_vec, _ = precon.apply(f_vec, image)
1345 precon_forces.append(pf_vec.reshape(-1, 3))
1347 return np.array(precon_forces)
1349 def average_norm(self, i, j, dx):
1350 """Average norm between images i and j
1352 Args:
1353 i (int): left image
1354 j (int): right image
1355 dx (array): vector
1357 Returns
1358 -------
1359 norm: norm of vector wrt average of precons at i and j
1360 """
1361 return np.sqrt(0.5 * (self.precon[i].dot(dx, dx) +
1362 self.precon[j].dot(dx, dx)))
1364 def get_tangent(self, i):
1365 """Normalised tangent vector at image i
1367 Args:
1368 i (int): image of interest
1370 Returns
1371 -------
1372 tangent: tangent vector, normalised with appropriate precon norm
1373 """
1374 tangent = self.spline.dx_ds(self.spline.s[i])
1375 tangent /= self.precon[i].norm(tangent)
1376 return tangent.reshape(-1, 3)
1378 def get_residual(self, i, imgforce):
1379 # residuals computed according to eq. 11 in the paper
1380 P_dot_imgforce = self.precon[i].Pdot(imgforce.reshape(-1))
1381 return np.linalg.norm(P_dot_imgforce, np.inf)
1383 def get_spring_force(self, i, k1, k2, tangent):
1384 """Spring force on image
1386 Args:
1387 i (int): image of interest
1388 k1 (float): spring constant for left spring
1389 k2 (float): spring constant for right spring
1390 tangent (array): tangent vector, shape (natoms, 3)
1392 Returns
1393 -------
1394 eta: NEB spring forces, shape (natoms, 3)
1395 """
1396 # Definition following Eq. 9 in Paper IV
1397 nimages = len(self.images)
1398 k = 0.5 * (k1 + k2) / (nimages ** 2)
1399 curvature = self.spline.d2x_ds2(self.spline.s[i]).reshape(-1, 3)
1400 # complete Eq. 9 by including the spring force
1401 eta = k * self.precon[i].vdot(curvature, tangent) * tangent
1402 return eta
1404 def get_coordinates(self, positions=None):
1405 """Compute displacements wrt appropriate precon metric for each image
1407 Args:
1408 positions (list or array, optional) - images positions.
1409 Shape either (nimages * natoms, 3) or ((nimages-2)*natoms, 3)
1411 Returns
1412 -------
1413 s : array shape (nimages,), reaction coordinates, in range [0, 1]
1414 x : array shape (nimages, 3 * natoms), flat displacement vectors
1415 """
1416 nimages = len(self.precon)
1417 natoms = len(self.images[0])
1418 d_P = np.zeros(nimages)
1419 x = np.zeros((nimages, 3 * natoms)) # flattened positions
1420 if positions is None:
1421 positions = [image.positions for image in self.images]
1422 elif isinstance(positions, np.ndarray) and len(positions.shape) == 2:
1423 positions = positions.reshape(-1, natoms, 3)
1424 positions = [positions[i, :, :] for i in range(len(positions))]
1425 if len(positions) == len(self.images) - 2:
1426 # prepend and append the non-moving images
1427 positions = ([self.images[0].positions] + positions +
1428 [self.images[-1].positions])
1429 assert len(positions) == len(self.images)
1431 x[0, :] = positions[0].reshape(-1)
1432 for i in range(1, nimages):
1433 x[i, :] = positions[i].reshape(-1)
1434 dx, _ = find_mic(positions[i] - positions[i - 1],
1435 self.images[i - 1].cell,
1436 self.images[i - 1].pbc)
1437 dx = dx.reshape(-1)
1438 d_P[i] = self.average_norm(i, i - 1, dx)
1440 s = d_P.cumsum() / d_P.sum() # Eq. A1 in paper IV
1441 return s, x
1443 def spline_fit(self, positions=None):
1444 """Fit 3 * natoms cubic splines as a function of reaction coordinate
1446 Returns
1447 -------
1448 fit : :class:`ase.optimize.precon.SplineFit` object
1449 """
1450 s, x = self.get_coordinates(positions)
1451 return SplineFit(s, x)
1453 @property
1454 def spline(self):
1455 s, x = self.get_coordinates()
1456 if self._spline and (np.abs(s - self._old_s).max() < 1e-6 and
1457 np.abs(x - self._old_x).max() < 1e-6):
1458 return self._spline
1460 self._spline = self.spline_fit()
1461 self._old_s = s
1462 self._old_x = x
1463 return self._spline