r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1# flake8: noqa

2"""Minimum mode follower for finding saddle points in an unbiased way.

4There is, currently, only one implemented method: The Dimer method.

6"""

8import sys

9import time

10import warnings

11from math import cos, sin, atan, tan, degrees, pi, sqrt

12from typing import Dict, Any

14import numpy as np

16from ase.optimize.optimize import Optimizer

17from ase.parallel import world

18from ase.calculators.singlepoint import SinglePointCalculator

19from ase.utils import IOContext

21# Handy vector methods

22norm = np.linalg.norm

25def normalize(vector):

26 """Create a unit vector along *vector*"""

27 return vector / norm(vector)

30def parallel_vector(vector, base):

31 """Extract the components of *vector* that are parallel to *base*"""

32 return np.vdot(vector, base) * base

35def perpendicular_vector(vector, base):

36 """Remove the components of *vector* that are parallel to *base*"""

37 return vector - parallel_vector(vector, base)

40def rotate_vectors(v1i, v2i, angle):

41 """Rotate vectors *v1i* and *v2i* by *angle*"""

42 cAng = cos(angle)

43 sAng = sin(angle)

44 v1o = v1i * cAng + v2i * sAng

45 v2o = v2i * cAng - v1i * sAng

46 # Ensure the length of the input and output vectors is equal

47 return normalize(v1o) * norm(v1i), normalize(v2o) * norm(v2i)

50class DimerEigenmodeSearch:

51 """An implementation of the Dimer's minimum eigenvalue mode search.

53 This class implements the rotational part of the dimer saddle point

54 searching method.

56 Parameters:

58 atoms: MinModeAtoms object

59 MinModeAtoms is an extension to the Atoms object, which includes

60 information about the lowest eigenvalue mode.

61 control: DimerControl object

62 Contains the parameters necessary for the eigenmode search.

63 If no control object is supplied a default DimerControl

64 will be created and used.

65 basis: list of xyz-values

66 Eigenmode. Must be an ndarray of shape (n, 3).

67 It is possible to constrain the eigenmodes to be orthogonal

68 to this given eigenmode.

70 Notes:

72 The code is inspired, with permission, by code written by the Henkelman

73 group, which can be found at http://theory.cm.utexas.edu/vtsttools/code/

75 References:

77 * Henkelman and Jonsson, JCP 111, 7010 (1999)

78 * Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,

79 9776 (2004).

80 * Heyden, Bell, and Keil, JCP 123, 224101 (2005).

81 * Kastner and Sherwood, JCP 128, 014106 (2008).

83 """

85 def __init__(self, atoms, control=None, eigenmode=None, basis=None,

86 **kwargs):

87 if hasattr(atoms, 'get_eigenmode'):

88 self.atoms = atoms

89 else:

90 e = 'The atoms object must be a MinModeAtoms object'

91 raise TypeError(e)

92 self.basis = basis

94 if eigenmode is None:

95 self.eigenmode = self.atoms.get_eigenmode()

96 else:

97 self.eigenmode = eigenmode

99 if control is None:

100 self.control = DimerControl(**kwargs)

101 w = 'Missing control object in ' + self.__class__.__name__ + \

102 '. Using default: DimerControl()'

103 warnings.warn(w, UserWarning)

104 if self.control.logfile is not None:

105 self.control.logfile.write('DIM:WARN: ' + w + '\n')

106 self.control.logfile.flush()

107 else:

108 self.control = control

109 # kwargs must be empty if a control object is supplied

110 for key in kwargs:

111 e = '__init__() got an unexpected keyword argument \'%s\'' % \

112 (key)

113 raise TypeError(e)

115 self.dR = self.control.get_parameter('dimer_separation')

116 self.logfile = self.control.get_logfile()

118 def converge_to_eigenmode(self):

119 """Perform an eigenmode search."""

120 self.set_up_for_eigenmode_search()

121 stoprot = False

123 # Load the relevant parameters from control

124 f_rot_min = self.control.get_parameter('f_rot_min')

125 f_rot_max = self.control.get_parameter('f_rot_max')

126 trial_angle = self.control.get_parameter('trial_angle')

127 max_num_rot = self.control.get_parameter('max_num_rot')

128 extrapolate = self.control.get_parameter('extrapolate_forces')

130 while not stoprot:

131 if self.forces1E is None:

132 self.update_virtual_forces()

133 else:

134 self.update_virtual_forces(extrapolated_forces=True)

135 self.forces1A = self.forces1

136 self.update_curvature()

137 f_rot_A = self.get_rotational_force()

139 # Pre rotation stop criteria

140 if norm(f_rot_A) <= f_rot_min:

141 self.log(f_rot_A, None)

142 stoprot = True

143 else:

144 n_A = self.eigenmode

145 rot_unit_A = normalize(f_rot_A)

147 # Get the curvature and its derivative

148 c0 = self.get_curvature()

149 c0d = np.vdot((self.forces2 - self.forces1), rot_unit_A) / \

150 self.dR

152 # Trial rotation (no need to store the curvature)

153 # NYI variable trial angles from [3]

154 n_B, rot_unit_B = rotate_vectors(n_A, rot_unit_A, trial_angle)

155 self.eigenmode = n_B

156 self.update_virtual_forces()

157 self.forces1B = self.forces1

159 # Get the curvature's derivative

160 c1d = np.vdot((self.forces2 - self.forces1), rot_unit_B) / \

161 self.dR

163 # Calculate the Fourier coefficients

164 a1 = c0d * cos(2 * trial_angle) - c1d / \

165 (2 * sin(2 * trial_angle))

166 b1 = 0.5 * c0d

167 a0 = 2 * (c0 - a1)

169 # Estimate the rotational angle

170 rotangle = atan(b1 / a1) / 2.0

172 # Make sure that you didn't find a maximum

173 cmin = a0 / 2.0 + a1 * cos(2 * rotangle) + \

174 b1 * sin(2 * rotangle)

175 if c0 < cmin:

176 rotangle += pi / 2.0

178 # Rotate into the (hopefully) lowest eigenmode

179 # NYI Conjugate gradient rotation

180 n_min, dummy = rotate_vectors(n_A, rot_unit_A, rotangle)

181 self.update_eigenmode(n_min)

183 # Store the curvature estimate instead of the old curvature

184 self.update_curvature(cmin)

186 self.log(f_rot_A, rotangle)

188 # Force extrapolation scheme from [4]

189 if extrapolate:

190 self.forces1E = sin(trial_angle - rotangle) / \

191 sin(trial_angle) * self.forces1A + sin(rotangle) / \

192 sin(trial_angle) * self.forces1B + \

193 (1 - cos(rotangle) - sin(rotangle) *

194 tan(trial_angle / 2.0)) * self.forces0

195 else:

196 self.forces1E = None

198 # Post rotation stop criteria

199 if not stoprot:

200 if self.control.get_counter('rotcount') >= max_num_rot:

201 stoprot = True

202 elif norm(f_rot_A) <= f_rot_max:

203 stoprot = True

205 def log(self, f_rot_A, angle):

206 """Log each rotational step."""

207 # NYI Log for the trial angle

208 if self.logfile is not None:

209 if angle:

210 l = 'DIM:ROT: %7d %9d %9.4f %9.4f %9.4f\n' % \

211 (self.control.get_counter('optcount'),

212 self.control.get_counter('rotcount'),

213 self.get_curvature(), degrees(angle), norm(f_rot_A))

214 else:

215 l = 'DIM:ROT: %7d %9d %9.4f %9s %9.4f\n' % \

216 (self.control.get_counter('optcount'),

217 self.control.get_counter('rotcount'),

218 self.get_curvature(), '---------', norm(f_rot_A))

219 self.logfile.write(l)

220 self.logfile.flush()

222 def get_rotational_force(self):

223 """Calculate the rotational force that acts on the dimer."""

224 rot_force = perpendicular_vector((self.forces1 - self.forces2),

225 self.eigenmode) / (2.0 * self.dR)

226 if self.basis is not None:

227 if len(self.basis) == len(self.atoms) and len(self.basis[0]) == \

228 3 and isinstance(self.basis[0][0], float):

229 rot_force = perpendicular_vector(rot_force, self.basis)

230 else:

231 for base in self.basis:

232 rot_force = perpendicular_vector(rot_force, base)

233 return rot_force

235 def update_curvature(self, curv=None):

236 """Update the curvature in the MinModeAtoms object."""

237 if curv:

238 self.curvature = curv

239 else:

240 self.curvature = np.vdot((self.forces2 - self.forces1),

241 self.eigenmode) / (2.0 * self.dR)

243 def update_eigenmode(self, eigenmode):

244 """Update the eigenmode in the MinModeAtoms object."""

245 self.eigenmode = eigenmode

246 self.update_virtual_positions()

247 self.control.increment_counter('rotcount')

249 def get_eigenmode(self):

250 """Returns the current eigenmode."""

251 return self.eigenmode

253 def get_curvature(self):

254 """Returns the curvature along the current eigenmode."""

255 return self.curvature

257 def get_control(self):

258 """Return the control object."""

259 return self.control

261 def update_center_forces(self):

262 """Get the forces at the center of the dimer."""

263 self.atoms.set_positions(self.pos0)

264 self.forces0 = self.atoms.get_forces(real=True)

265 self.energy0 = self.atoms.get_potential_energy()

267 def update_virtual_forces(self, extrapolated_forces=False):

268 """Get the forces at the endpoints of the dimer."""

269 self.update_virtual_positions()

271 # Estimate / Calculate the forces at pos1

272 if extrapolated_forces:

273 self.forces1 = self.forces1E.copy()

274 else:

275 self.forces1 = self.atoms.get_forces(real=True, pos=self.pos1)

277 # Estimate / Calculate the forces at pos2

278 if self.control.get_parameter('use_central_forces'):

279 self.forces2 = 2 * self.forces0 - self.forces1

280 else:

281 self.forces2 = self.atoms.get_forces(real=True, pos=self.pos2)

283 def update_virtual_positions(self):

284 """Update the end point positions."""

285 self.pos1 = self.pos0 + self.eigenmode * self.dR

286 self.pos2 = self.pos0 - self.eigenmode * self.dR

288 def set_up_for_eigenmode_search(self):

289 """Before eigenmode search, prepare for rotation."""

290 self.pos0 = self.atoms.get_positions()

291 self.update_center_forces()

292 self.update_virtual_positions()

293 self.control.reset_counter('rotcount')

294 self.forces1E = None

296 def set_up_for_optimization_step(self):

297 """At the end of rotation, prepare for displacement of the dimer."""

298 self.atoms.set_positions(self.pos0)

299 self.forces1E = None

302class MinModeControl(IOContext):

303 """A parent class for controlling minimum mode saddle point searches.

305 Method specific control classes inherit this one. The only thing

306 inheriting classes need to implement are the log() method and

307 the *parameters* class variable with default values for ALL

308 parameters needed by the method in question.

309 When instantiating control classes default parameter values can

310 be overwritten.

312 """

313 parameters: Dict[str, Any] = {}

315 def __init__(self, logfile='-', eigenmode_logfile=None, **kwargs):

316 # Overwrite the defaults with the input parameters given

317 for key in kwargs:

318 if not key in self.parameters:

319 e = 'Invalid parameter >>%s<< with value >>%s<< in %s' % \

320 (key, str(kwargs[key]), self.__class__.__name__)

321 raise ValueError(e)

322 else:

323 self.set_parameter(key, kwargs[key], log=False)

325 self.initialize_logfiles(logfile, eigenmode_logfile)

326 self.counters = {'forcecalls': 0, 'rotcount': 0, 'optcount': 0}

327 self.log()

329 def initialize_logfiles(self, logfile=None, eigenmode_logfile=None):

330 self.logfile = self.openfile(logfile, comm=world)

331 self.eigenmode_logfile = self.openfile(eigenmode_logfile, comm=world)

333 def log(self, parameter=None):

334 """Log the parameters of the eigenmode search."""

335 pass

337 def set_parameter(self, parameter, value, log=True):

338 """Change a parameter's value."""

339 if not parameter in self.parameters:

340 e = 'Invalid parameter >>%s<< with value >>%s<<' % \

341 (parameter, str(value))

342 raise ValueError(e)

343 self.parameters[parameter] = value

344 if log:

345 self.log(parameter)

347 def get_parameter(self, parameter):

348 """Returns the value of a parameter."""

349 if not parameter in self.parameters:

350 e = 'Invalid parameter >>%s<<' % \

351 (parameter)

352 raise ValueError(e)

353 return self.parameters[parameter]

355 def get_logfile(self):

356 """Returns the log file."""

357 return self.logfile

359 def get_eigenmode_logfile(self):

360 """Returns the eigenmode log file."""

361 return self.eigenmode_logfile

363 def get_counter(self, counter):

364 """Returns a given counter."""

365 return self.counters[counter]

367 def increment_counter(self, counter):

368 """Increment a given counter."""

369 self.counters[counter] += 1

371 def reset_counter(self, counter):

372 """Reset a given counter."""

373 self.counters[counter] = 0

375 def reset_all_counters(self):

376 """Reset all counters."""

377 for key in self.counters:

378 self.counters[key] = 0

381class DimerControl(MinModeControl):

382 """A class that takes care of the parameters needed for a Dimer search.

384 Parameters:

386 eigenmode_method: str

387 The name of the eigenmode search method.

388 f_rot_min: float

389 Size of the rotational force under which no rotation will be

390 performed.

391 f_rot_max: float

392 Size of the rotational force under which only one rotation will be

393 performed.

394 max_num_rot: int

395 Maximum number of rotations per optimizer step.

396 trial_angle: float

397 Trial angle for the finite difference estimate of the rotational

399 trial_trans_step: float

400 Trial step size for the MinModeTranslate optimizer.

401 maximum_translation: float

402 Maximum step size and forced step size when the curvature is still

403 positive for the MinModeTranslate optimizer.

404 cg_translation: bool

405 Conjugate Gradient for the MinModeTranslate optimizer.

406 use_central_forces: bool

407 Only calculate the forces at one end of the dimer and extrapolate

408 the forces to the other.

409 dimer_separation: float

410 Separation of the dimer's images.

411 initial_eigenmode_method: str

412 How to construct the initial eigenmode of the dimer. If an eigenmode

413 is given when creating the MinModeAtoms object, this will be ignored.

414 Possible choices are: 'gauss' and 'displacement'

415 extrapolate_forces: bool

416 When more than one rotation is performed, an extrapolation scheme can

417 be used to reduce the number of force evaluations.

418 displacement_method: str

419 How to displace the atoms. Possible choices are 'gauss' and 'vector'.

420 gauss_std: float

421 The standard deviation of the gauss curve used when doing random

422 displacement.

423 order: int

424 How many lowest eigenmodes will be inverted.

426 Which atoms will be moved during displacement.

427 displacement_center: int or [float, float, float]

428 The center of displacement, nearby atoms will be displaced.

430 When choosing which atoms to displace with the *displacement_center*

431 keyword, this decides how many nearby atoms to displace.

432 number_of_displacement_atoms: int

433 The amount of atoms near *displacement_center* to displace.

435 """

436 # Default parameters for the Dimer eigenmode search

437 parameters = {'eigenmode_method': 'dimer',

438 'f_rot_min': 0.1,

439 'f_rot_max': 1.00,

440 'max_num_rot': 1,

441 'trial_angle': pi / 4.0,

442 'trial_trans_step': 0.001,

443 'maximum_translation': 0.1,

444 'cg_translation': True,

445 'use_central_forces': True,

446 'dimer_separation': 0.0001,

447 'initial_eigenmode_method': 'gauss',

448 'extrapolate_forces': False,

449 'displacement_method': 'gauss',

450 'gauss_std': 0.1,

451 'order': 1,

453 'displacement_center': None,

455 'number_of_displacement_atoms': None}

457 # NB: Can maybe put this in EigenmodeSearch and MinModeControl

458 def log(self, parameter=None):

459 """Log the parameters of the eigenmode search."""

460 if self.logfile is not None:

461 if parameter is not None:

462 l = 'DIM:CONTROL: Updated Parameter: %s = %s\n' % (parameter,

463 str(self.get_parameter(parameter)))

464 else:

465 l = 'MINMODE:METHOD: Dimer\n'

466 l += 'DIM:CONTROL: Search Parameters:\n'

467 l += 'DIM:CONTROL: ------------------\n'

468 for key in self.parameters:

469 l += 'DIM:CONTROL: %s = %s\n' % (key,

470 str(self.get_parameter(key)))

471 l += 'DIM:CONTROL: ------------------\n'

472 l += 'DIM:ROT: OPT-STEP ROT-STEP CURVATURE ROT-ANGLE ' + \

473 'ROT-FORCE\n'

474 self.logfile.write(l)

475 self.logfile.flush()

478class MinModeAtoms:

479 """Wrapper for Atoms with information related to minimum mode searching.

481 Contains an Atoms object and pipes all unknown function calls to that

482 object.

483 Other information that is stored in this object are the estimate for

484 the lowest eigenvalue, *curvature*, and its corresponding eigenmode,

485 *eigenmode*. Furthermore, the original configuration of the Atoms

486 object is stored for use in multiple minimum mode searches.

487 The forces on the system are modified by inverting the component

488 along the eigenmode estimate. This eventually brings the system to

491 Parameters:

493 atoms : Atoms object

494 A regular Atoms object

495 control : MinModeControl object

496 Contains the parameters necessary for the eigenmode search.

497 If no control object is supplied a default DimerControl

498 will be created and used.

500 Determines which atoms will be moved when calling displace()

501 random_seed: int

502 The seed used for the random number generator. Defaults to

503 modified version the current time.

505 References: [1]_ [2]_ [3]_ [4]_

507 .. [1] Henkelman and Jonsson, JCP 111, 7010 (1999)

508 .. [2] Olsen, Kroes, Henkelman, Arnaldsson, and Jonsson, JCP 121,

509 9776 (2004).

510 .. [3] Heyden, Bell, and Keil, JCP 123, 224101 (2005).

511 .. [4] Kastner and Sherwood, JCP 128, 014106 (2008).

513 """

515 def __init__(self, atoms, control=None, eigenmodes=None,

516 random_seed=None, **kwargs):

517 self.minmode_init = True

518 self.atoms = atoms

520 # Initialize to None to avoid strange behaviour due to __getattr__

521 self.eigenmodes = eigenmodes

522 self.curvatures = None

524 if control is None:

525 self.control = DimerControl(**kwargs)

526 w = 'Missing control object in ' + self.__class__.__name__ + \

527 '. Using default: DimerControl()'

528 warnings.warn(w, UserWarning)

529 if self.control.logfile is not None:

530 self.control.logfile.write('DIM:WARN: ' + w + '\n')

531 self.control.logfile.flush()

532 else:

533 self.control = control

534 logfile = self.control.get_logfile()

535 mlogfile = self.control.get_eigenmode_logfile()

536 for key in kwargs:

537 if key == 'logfile':

538 logfile = kwargs[key]

539 elif key == 'eigenmode_logfile':

540 mlogfile = kwargs[key]

541 else:

542 self.control.set_parameter(key, kwargs[key])

543 self.control.initialize_logfiles(logfile=logfile,

544 eigenmode_logfile=mlogfile)

546 # Seed the randomness

547 if random_seed is None:

548 t = time.time()

549 if world.size > 1:

550 t = world.sum(t) / world.size

551 # Harvest the latter part of the current time

552 random_seed = int(('%30.9f' % t)[-9:])

553 self.random_state = np.random.RandomState(random_seed)

555 # Check the order

556 self.order = self.control.get_parameter('order')

558 # Construct the curvatures list

559 self.curvatures = [100.0] * self.order

561 # Save the original state of the atoms.

562 self.atoms0 = self.atoms.copy()

563 self.save_original_forces()

565 # Get a reference to the log files

566 self.logfile = self.control.get_logfile()

567 self.mlogfile = self.control.get_eigenmode_logfile()

569 def save_original_forces(self, force_calculation=False):

570 """Store the forces (and energy) of the original state."""

571 # NB: Would be nice if atoms.copy() took care of this.

572 if self.calc is not None:

573 # Hack because some calculators do not have calculation_required

574 if (hasattr(self.calc, 'calculation_required')

575 and not self.calc.calculation_required(self.atoms,

576 ['energy', 'forces'])) or force_calculation:

577 calc = SinglePointCalculator(

578 self.atoms0,

579 energy=self.atoms.get_potential_energy(),

580 forces=self.atoms.get_forces())

581 self.atoms0.calc = calc

583 def initialize_eigenmodes(self, method=None, eigenmodes=None,

584 gauss_std=None):

585 """Make an initial guess for the eigenmode."""

586 if eigenmodes is None:

587 pos = self.get_positions()

588 old_pos = self.get_original_positions()

589 if method is None:

590 method = \

591 self.control.get_parameter('initial_eigenmode_method')

592 if method.lower() == 'displacement' and (pos - old_pos).any():

593 eigenmode = normalize(pos - old_pos)

594 elif method.lower() == 'gauss':

595 self.displace(log=False, gauss_std=gauss_std,

596 method=method)

597 new_pos = self.get_positions()

598 eigenmode = normalize(new_pos - pos)

599 self.set_positions(pos)

600 else:

601 e = 'initial_eigenmode must use either \'gauss\' or ' + \

602 '\'displacement\', if the latter is used the atoms ' + \

603 'must have moved away from the original positions.' + \

604 'You have requested \'%s\'.' % method

605 raise NotImplementedError(e) # NYI

606 eigenmodes = [eigenmode]

608 # Create random higher order mode guesses

609 if self.order > 1:

610 if len(eigenmodes) == 1:

611 for k in range(1, self.order):

612 pos = self.get_positions()

613 self.displace(log=False, gauss_std=gauss_std,

614 method=method)

615 new_pos = self.get_positions()

616 eigenmode = normalize(new_pos - pos)

617 self.set_positions(pos)

618 eigenmodes += [eigenmode]

620 self.eigenmodes = eigenmodes

621 # Ensure that the higher order mode guesses are all orthogonal

622 if self.order > 1:

623 for k in range(self.order):

624 self.ensure_eigenmode_orthogonality(k)

625 self.eigenmode_log()

627 # NB maybe this name might be confusing in context to

628 # calc.calculation_required()

629 def calculation_required(self):

630 """Check if a calculation is required."""

631 return self.minmode_init or self.check_atoms != self.atoms

633 def calculate_real_forces_and_energies(self, **kwargs):

634 """Calculate and store the potential energy and forces."""

635 if self.minmode_init:

636 self.minmode_init = False

637 self.initialize_eigenmodes(eigenmodes=self.eigenmodes)

638 self.rotation_required = True

639 self.forces0 = self.atoms.get_forces(**kwargs)

640 self.energy0 = self.atoms.get_potential_energy()

641 self.control.increment_counter('forcecalls')

642 self.check_atoms = self.atoms.copy()

644 def get_potential_energy(self):

645 """Return the potential energy."""

646 if self.calculation_required():

647 self.calculate_real_forces_and_energies()

648 return self.energy0

650 def get_forces(self, real=False, pos=None, **kwargs):

651 """Return the forces, projected or real."""

652 if self.calculation_required() and pos is None:

653 self.calculate_real_forces_and_energies(**kwargs)

654 if real and pos is None:

655 return self.forces0

656 elif real and pos is not None:

657 old_pos = self.atoms.get_positions()

658 self.atoms.set_positions(pos)

659 forces = self.atoms.get_forces()

660 self.control.increment_counter('forcecalls')

661 self.atoms.set_positions(old_pos)

662 return forces

663 else:

664 if self.rotation_required:

665 self.find_eigenmodes(order=self.order)

666 self.eigenmode_log()

667 self.rotation_required = False

668 self.control.increment_counter('optcount')

669 return self.get_projected_forces()

671 def ensure_eigenmode_orthogonality(self, order):

672 mode = self.eigenmodes[order - 1].copy()

673 for k in range(order - 1):

674 mode = perpendicular_vector(mode, self.eigenmodes[k])

675 self.eigenmodes[order - 1] = normalize(mode)

677 def find_eigenmodes(self, order=1):

678 """Launch eigenmode searches."""

679 if self.control.get_parameter('eigenmode_method').lower() != 'dimer':

680 e = 'Only the Dimer control object has been implemented.'

681 raise NotImplementedError(e) # NYI

682 for k in range(order):

683 if k > 0:

684 self.ensure_eigenmode_orthogonality(k + 1)

685 search = DimerEigenmodeSearch(self, self.control,

686 eigenmode=self.eigenmodes[k], basis=self.eigenmodes[:k])

687 search.converge_to_eigenmode()

688 search.set_up_for_optimization_step()

689 self.eigenmodes[k] = search.get_eigenmode()

690 self.curvatures[k] = search.get_curvature()

692 def get_projected_forces(self, pos=None):

693 """Return the projected forces."""

694 if pos is not None:

695 forces = self.get_forces(real=True, pos=pos).copy()

696 else:

697 forces = self.forces0.copy()

699 # Loop through all the eigenmodes

700 # NB: Can this be done with a linear combination, instead?

701 for k, mode in enumerate(self.eigenmodes):

702 # NYI This If statement needs to be overridable in the control

703 if self.get_curvature(order=k) > 0.0 and self.order == 1:

704 forces = -parallel_vector(forces, mode)

705 else:

706 forces -= 2 * parallel_vector(forces, mode)

707 return forces

709 def restore_original_positions(self):

710 """Restore the MinModeAtoms object positions to the original state."""

711 self.atoms.set_positions(self.get_original_positions())

713 def get_barrier_energy(self):

714 """The energy difference between the current and original states"""

715 try:

716 original_energy = self.get_original_potential_energy()

717 dimer_energy = self.get_potential_energy()

718 return dimer_energy - original_energy

719 except RuntimeError:

720 w = 'The potential energy is not available, without further ' + \

721 'calculations, most likely at the original state.'

722 warnings.warn(w, UserWarning)

723 return np.nan

725 def get_control(self):

726 """Return the control object."""

727 return self.control

729 def get_curvature(self, order='max'):

730 """Return the eigenvalue estimate."""

731 if order == 'max':

732 return max(self.curvatures)

733 else:

734 return self.curvatures[order - 1]

736 def get_eigenmode(self, order=1):

737 """Return the current eigenmode guess."""

738 return self.eigenmodes[order - 1]

740 def get_atoms(self):

741 """Return the unextended Atoms object."""

742 return self.atoms

744 def set_atoms(self, atoms):

745 """Set a new Atoms object"""

746 self.atoms = atoms

748 def set_eigenmode(self, eigenmode, order=1):

749 """Set the eigenmode guess."""

750 self.eigenmodes[order - 1] = eigenmode

752 def set_curvature(self, curvature, order=1):

753 """Set the eigenvalue estimate."""

754 self.curvatures[order - 1] = curvature

756 # Pipe all the stuff from Atoms that is not overwritten.

757 # Pipe all requests for get_original_* to self.atoms0.

758 def __getattr__(self, attr):

759 """Return any value of the Atoms object"""

760 if 'original' in attr.split('_'):

761 attr = attr.replace('_original_', '_')

762 return getattr(self.atoms0, attr)

763 else:

764 return getattr(self.atoms, attr)

766 def __len__(self):

767 return len(self.atoms)

769 def displace(self, displacement_vector=None, mask=None, method=None,

771 gauss_std=None, mic=True, log=True):

772 """Move the atoms away from their current position.

774 This is one of the essential parts of minimum mode searches.

775 The parameters can all be set in the control object and overwritten

776 when this method is run, apart from *displacement_vector*.

777 It is preferred to modify the control values rather than those here

778 in order for the correct ones to show up in the log file.

780 *method* can be either 'gauss' for random displacement or 'vector'

781 to perform a predefined displacement.

783 *gauss_std* is the standard deviation of the gauss curve that is

784 used for random displacement.

786 *displacement_center* can be either the number of an atom or a 3D

787 position. It must be accompanied by a *radius* (all atoms within it

788 will be displaced) or a *number_of_atoms* which decides how many of

789 the closest atoms will be displaced.

791 *mic* controls the usage of the Minimum Image Convention.

793 If both *mask* and *displacement_center* are used, the atoms marked

794 as False in the *mask* will not be affected even though they are

795 within reach of the *displacement_center*.

797 The parameters priority order:

798 1) displacement_vector

800 3) displacement_center (with radius and/or number_of_atoms)

802 If both *radius* and *number_of_atoms* are supplied with

803 *displacement_center*, only atoms that fulfill both criteria will

804 be displaced.

806 """

808 # Fetch the default values from the control

811 if method is None:

812 method = self.control.get_parameter('displacement_method')

813 if gauss_std is None:

814 gauss_std = self.control.get_parameter('gauss_std')

815 if displacement_center is None:

816 displacement_center = \

817 self.control.get_parameter('displacement_center')

820 if number_of_atoms is None:

821 number_of_atoms = \

822 self.control.get_parameter('number_of_displacement_atoms')

824 # Check for conflicts

825 if displacement_vector is not None and method.lower() != 'vector':

826 e = 'displacement_vector was supplied but a different method ' + \

827 '(\'%s\') was chosen.\n' % str(method)

828 raise ValueError(e)

829 elif displacement_vector is None and method.lower() == 'vector':

830 e = 'A displacement_vector must be supplied when using ' + \

831 'method = \'%s\'.\n' % str(method)

832 raise ValueError(e)

833 elif displacement_center is not None and radius is None and \

834 number_of_atoms is None:

835 e = 'When displacement_center is chosen, either radius or ' + \

836 'number_of_atoms must be supplied.\n'

837 raise ValueError(e)

840 if displacement_center is not None:

841 c = displacement_center

842 # Construct a distance list

843 # The center is an atom

844 if isinstance(c, int):

845 # Parse negative indexes

846 c = displacement_center % len(self)

847 d = [(k, self.get_distance(k, c, mic=mic)) for k in

848 range(len(self))]

849 # The center is a position in 3D space

850 elif len(c) == 3 and [type(c_k) for c_k in c] == [float] * 3:

851 # NB: MIC is not considered.

852 d = [(k, norm(self.get_positions()[k] - c))

853 for k in range(len(self))]

854 else:

855 e = 'displacement_center must be either the number of an ' + \

856 'atom in MinModeAtoms object or a 3D position ' + \

857 '(3-tuple of floats).'

858 raise ValueError(e)

860 # Set up the mask

861 if radius is not None:

863 else:

864 r_mask = [True for _ in range(len(self))]

866 if number_of_atoms is not None:

867 d_sorted = [n[0] for n in sorted(d, key=lambda k: k[1])]

868 n_nearest = d_sorted[:number_of_atoms]

869 n_mask = [k in n_nearest for k in range(len(self))]

870 else:

871 n_mask = [True for _ in range(len(self))]

875 else:

878 # Set up a True mask if there is no mask supplied

880 mask = [True for _ in range(len(self))]

882 w = 'It was not possible to figure out which atoms to ' + \

883 'displace, Will try to displace all atoms.\n'

884 warnings.warn(w, UserWarning)

885 if self.logfile is not None:

886 self.logfile.write('MINMODE:WARN: ' + w + '\n')

887 self.logfile.flush()

890 if c_mask is not None:

893 if displacement_vector is None:

894 displacement_vector = []

895 for k in range(len(self)):

897 diff_line = []

898 for _ in range(3):

899 if method.lower() == 'gauss':

900 if not gauss_std:

901 gauss_std = \

902 self.control.get_parameter('gauss_std')

903 diff = self.random_state.normal(0.0, gauss_std)

904 else:

905 e = 'Invalid displacement method >>%s<<' % \

906 str(method)

907 raise ValueError(e)

908 diff_line.append(diff)

909 displacement_vector.append(diff_line)

910 else:

911 displacement_vector.append([0.0] * 3)

913 # Remove displacement of masked atoms

916 displacement_vector[k] = [0.0] * 3

918 # Perform the displacement and log it

919 if log:

920 pos0 = self.get_positions()

921 self.set_positions(self.get_positions() + displacement_vector)

922 if log:

924 'displacement_method': method,

925 'gauss_std': gauss_std,

926 'displacement_center': displacement_center,

928 'number_of_displacement_atoms': number_of_atoms}

929 self.displacement_log(self.get_positions() - pos0, parameters)

931 def eigenmode_log(self):

932 """Log the eigenmodes (eigenmode estimates)"""

933 if self.mlogfile is not None:

934 l = 'MINMODE:MODE: Optimization Step: %i\n' % \

935 (self.control.get_counter('optcount'))

936 for m_num, mode in enumerate(self.eigenmodes):

937 l += 'MINMODE:MODE: Order: %i\n' % m_num

938 for k in range(len(mode)):

939 l += 'MINMODE:MODE: %7i %15.8f %15.8f %15.8f\n' % (k,

940 mode[k][0], mode[k][1], mode[k][2])

941 self.mlogfile.write(l)

942 self.mlogfile.flush()

944 def displacement_log(self, displacement_vector, parameters):

945 """Log the displacement"""

946 if self.logfile is not None:

947 lp = 'MINMODE:DISP: Parameters, different from the control:\n'

948 mod_para = False

949 for key in parameters:

950 if parameters[key] != self.control.get_parameter(key):

951 lp += 'MINMODE:DISP: %s = %s\n' % (str(key),

952 str(parameters[key]))

953 mod_para = True

954 if mod_para:

955 l = lp

956 else:

957 l = ''

958 for k in range(len(displacement_vector)):

959 l += 'MINMODE:DISP: %7i %15.8f %15.8f %15.8f\n' % (k,

960 displacement_vector[k][0], displacement_vector[k][1],

961 displacement_vector[k][2])

962 self.logfile.write(l)

963 self.logfile.flush()

965 def summarize(self):

966 """Summarize the Minimum mode search."""

967 if self.logfile is None:

968 logfile = sys.stdout

969 else:

970 logfile = self.logfile

972 c = self.control

973 label = 'MINMODE:SUMMARY: '

975 l = label + '-------------------------\n'

976 l += label + 'Barrier: %16.4f\n' % self.get_barrier_energy()

977 l += label + 'Curvature: %14.4f\n' % self.get_curvature()

978 l += label + 'Optimizer steps: %8i\n' % c.get_counter('optcount')

979 l += label + 'Forcecalls: %13i\n' % c.get_counter('forcecalls')

980 l += label + '-------------------------\n'

982 logfile.write(l)

985class MinModeTranslate(Optimizer):

986 """An Optimizer specifically tailored to minimum mode following."""

988 def __init__(self, atoms, logfile='-', trajectory=None):

989 Optimizer.__init__(self, atoms, None, logfile, trajectory)

991 self.control = atoms.get_control()

993 # Make a header for the log

994 if self.logfile is not None:

995 l = ''

996 if isinstance(self.control, DimerControl):

997 l = 'MinModeTranslate: STEP TIME ENERGY ' + \

998 'MAX-FORCE STEPSIZE CURVATURE ROT-STEPS\n'

999 self.logfile.write(l)

1000 self.logfile.flush()

1002 # Load the relevant parameters from control

1003 self.cg_on = self.control.get_parameter('cg_translation')

1004 self.trial_step = self.control.get_parameter('trial_trans_step')

1005 self.max_step = self.control.get_parameter('maximum_translation')

1008 if self.cg_on:

1009 self.cg_init = True

1011 def initialize(self):

1012 """Set initial values."""

1013 self.r0 = None

1014 self.f0 = None

1016 def step(self, f=None):

1017 """Perform the optimization step."""

1018 atoms = self.atoms

1019 if f is None:

1020 f = atoms.get_forces()

1021 r = atoms.get_positions()

1022 curv = atoms.get_curvature()

1023 f0p = f.copy()

1024 r0 = r.copy()

1025 direction = f0p.copy()

1026 if self.cg_on:

1027 direction = self.get_cg_direction(direction)

1028 direction = normalize(direction)

1029 if curv > 0.0:

1030 step = direction * self.max_step

1031 else:

1032 r0t = r0 + direction * self.trial_step

1033 f0tp = self.atoms.get_projected_forces(r0t)

1034 F = np.vdot((f0tp + f0p), direction) / 2.0

1035 C = np.vdot((f0tp - f0p), direction) / self.trial_step

1036 step = (-F / C + self.trial_step / 2.0) * direction

1037 if norm(step) > self.max_step:

1038 step = direction * self.max_step

1039 self.log(f0p, norm(step))

1041 atoms.set_positions(r + step)

1043 self.f0 = f.flat.copy()

1044 self.r0 = r.flat.copy()

1046 def get_cg_direction(self, direction):

1047 """Apply the Conjugate Gradient algorithm to the step direction."""

1048 if self.cg_init:

1049 self.cg_init = False

1050 self.direction_old = direction.copy()

1051 self.cg_direction = direction.copy()

1052 old_norm = np.vdot(self.direction_old, self.direction_old)

1054 if old_norm != 0.0:

1055 betaPR = np.vdot(direction, (direction - self.direction_old)) / \

1056 old_norm

1057 else:

1058 betaPR = 0.0

1059 if betaPR < 0.0:

1060 betaPR = 0.0

1061 self.cg_direction = direction + self.cg_direction * betaPR

1062 self.direction_old = direction.copy()

1063 return self.cg_direction.copy()

1065 def log(self, f=None, stepsize=None):

1066 """Log each step of the optimization."""

1067 if f is None:

1068 f = self.atoms.get_forces()

1069 if self.logfile is not None:

1070 T = time.localtime()

1071 e = self.atoms.get_potential_energy()

1072 fmax = sqrt((f**2).sum(axis=1).max())

1073 rotsteps = self.atoms.control.get_counter('rotcount')

1074 curvature = self.atoms.get_curvature()

1075 l = ''

1076 if stepsize:

1077 if isinstance(self.control, DimerControl):

1078 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %12.6f ' \

1079 '%12.6f %10d\n' % ('MinModeTranslate', self.nsteps,

1080 T[3], T[4], T[5], e, fmax, stepsize, curvature,

1081 rotsteps)

1082 else:

1083 if isinstance(self.control, DimerControl):

1084 l = '%s: %4d %02d:%02d:%02d %15.6f %12.4f %s ' \

1085 '%12.6f %10d\n' % ('MinModeTranslate', self.nsteps,

1086 T[3], T[4], T[5], e, fmax, ' --------',

1087 curvature, rotsteps)

1088 self.logfile.write(l)

1089 self.logfile.flush()

1094 To access the pre optimization eigenmode set index = 'null'.

1096 """

1097 mlog_is_str = isinstance(mlog, str)

1098 if mlog_is_str:

1099 fd = open(mlog, 'r')

1100 else:

1101 fd = mlog

1105 # Detect the amount of atoms and iterations

1106 k = 2

1107 while lines[k].split()[1].lower() not in ['optimization', 'order']:

1108 k += 1

1109 n = k - 2

1110 n_itr = (len(lines) // (n + 1)) - 2

1112 # Locate the correct image.

1113 if isinstance(index, str):

1114 if index.lower() == 'null':

1115 i = 0

1116 else:

1117 i = int(index) + 1

1118 else:

1119 if index >= 0:

1120 i = index + 1

1121 else:

1122 if index < -n_itr - 1:

1123 raise IndexError('list index out of range')

1124 else:

1125 i = index

1127 mode = np.ndarray(shape=(n, 3), dtype=float)

1128 k_atom = 0

1129 for k in range(1, n + 1):

1130 line = lines[i * (n + 1) + k].split()

1131 for k_dim in range(3):

1132 mode[k_atom][k_dim] = float(line[k_dim + 2])

1133 k_atom += 1

1135 if mlog_is_str:

1136 fd.close()

1138 return mode

1141# Aliases

1142DimerAtoms = MinModeAtoms

1143DimerTranslate = MinModeTranslate