Algorithm¶
The optimization algorithm implemented in PyBerny loosely follows the “standard method” (SM) described in the appendix of [BirkholzTCA16]. What follows is a summary of that method, more detailed specification when necessary, and description of any deviations. For a full reference-manual treatment of the SM itself, with every equation needed to reimplement it and an explicit catalogue of PyBerny’s deviations, see Standard method — full reference.
Todo
Make the algorithm fully conform to the SM.
Sketch of the algorithm¶
Form a diagonal Hessian guess [SwartIJQC06].
Obtain energy and Cartesian gradients for the current geometry.
Form the Wilson B matrix and its generalized inverse.
Update the Hessian using the BFGS scheme.
Update trust region (Eq. 5.1.6 of [Fletcher00]).
Perform linear search (Gaussian website, section Examples, “If a minimum is sought…”). (Not in the SM.)
Project to a nonredundant subspace [PengJCC96].
Perform a quadratic RFO step [BanerjeeJPC85].
Transform back to Cartesian coordinates [PengJCC96].
If convergence is not reached (criteria from the SM), go to 3.
Symmetric start geometries¶
A gradient-following optimizer cannot leave the symmetric subspace of an exactly
symmetric start geometry: the gradient along every non-totally-symmetric mode
vanishes by symmetry, so the optimization can converge to a symmetric saddle
rather than a minimum (issue #148). By default Berny detects the
point group of the start and warns when it is not C1; pass
symmetry='break' to displace the start off its symmetry elements with a
small, deterministic, symmetry-targeted kick (the equal-weight sum of the
non-totally-symmetric Cartesian symmetry-adapted displacement coordinates) so
the optimizer can relax to the true minimum.
Redundant internal coordinates¶
All bonds shorter than 1.3 times the sum of covalent radii are created [PengJCC96].
If there are unconnected fragments, all bonds between unconnected fragments shorter than the sum of van der Waals radii plus d are created, with d starting at 0 and increasing by 1 angstrom, until all fragments are connected (custom scheme by JH).
All angles greater than 45° are created.
All dihedrals with 1–2–3, 2–3–4 angles both greater than 45° are created. If one of the angles is zero, so that three atoms lie on a line, they are used as a new base for a dihedral. This process is recursively repeated [PengJCC96].
In the case of a crystal, just the internal coordinate closest to the original unit cell is retained from all its periodic images.
Todo
Implement linear bends.
Generalized inverse¶
The Wilson B matrix, which relates differences in the internal redundant coordinates to differences in the Cartesian coordinates, is in general non-square and non-invertible. Its generalized inverse is obtained from the pseudoinverse of (singular), which is in turn obtained via singular value decomposition and inversion of only the nonzero singular values. For invertible matrices, this procedure is equivalent to an ordinary inverse. In practice, the zero values are in fact nonzero but several orders of magnitude smaller than the true nonzero values.
References¶
Birkholz, A. B. & Schlegel, H. B. Exploration of some refinements to geometry optimization methods. Theor. Chem. Acc. 135, (2016). DOI: 10.1007/s00214-016-1847-3
Peng, C., Ayala, P. Y., Schlegel, H. B. & Frisch, M. J. Using redundant internal coordinates to optimize equilibrium geometries and transition states. J. Comput. Chem. 17, 49–56 (1996). DOI: 10.1002/(SICI)1096-987X(19960115)17:1<49::AID-JCC5>3.0.CO;2-0
Swart, M. & Bickelhaupt, F. M. Optimization of Strong and Weak Coordinates. Int. J. Quantum Chem. 106, 2536–2544 (2006). DOI: 10.1002/qua.21049
Fletcher, R. Practical Methods of Optimization. (Wiley, 2000). URL: https://www.wiley.com/en-us/Practical+Methods+of+Optimization%2C+2nd+Edition-p-9780471494638
Banerjee, A., Adams, N. & Simons, J. Search for Stationary Points on Surfaces. J. Phys. Chem. 57, 52–57 (1985). DOI: 10.1021/j100247a015