8. Minimization with Microiterations
Pysisyphus allows layered optimizations using microiterations via the LayerOpt class. Required energies and its derivatives (gradients, Hessians) can either be calculated with pysisyphus' own ONIOM implementation or supplied through Unix sockets via the i-PI protocol.
In ONIOM calculations a small (model) system is treated using a high level of theory, while its surroundings (real system) are described at a lower, more economic level of theory. There may also be several (optional) intermediate layers. Please see the figure below.
Fig. 8.1 General structure of a n-layer ONIOM calculation in pysisyphus. Every box corresponds to a layer. Indices of the respective layers are given as red numbers in the top right corner of each layer.
Searching stationary points of such systems using classical optimization approaches with simultanous relaxation of model and real system can be computationally quite costly. In cases where the model system converges before the real system, many unnecessary model system energy and gradients evaluations may be required to fully relax the total system, as one step requires gradient evaluations in all layers. Also, it may be prohibitively expensive to optimize the whole system using internal coordinates, so it would be desireable to treat only the model system using internal coordinates.
Optimization with microiterations offers a way to potentially reduce the number of costly calculations in the innermost layer. By fully relaxing the outer layers before taking a step in the innermost layer, superfluous gradient evaluations can be avoided. As the overall calculation time of an ONIOM energy/gradient is usually dominated by the innermost layer, great compuational saving can by realized by decreasing the required number of such calculations.
In the present approach, outer layers are optimized in Cartesian coordinates with an economic optimizer like (preconditioned) limited memory BFGS (L-BFGS) or conjugate gradient (CG). As L-BFGS avoids operations with steep scaling costs like matrix diagonalization, relaxation of large system comprising thousands of atoms becomes possible with negligible computational overhead. In contrast to the outer layers, the inner layer is usually optimized using internal coordinates and an optimizer that utilizes an explicit Hessian matrix, as this often allows for the most efficient optimizations.
As energy and gradient evaluations in the outer layers should be cheap the inferior performance of optimizations in Cartesian coordinates should not lead to an overall runtime increase.
8.1. YAML input
When pysisyphus' native ONIOM calculator is to be used, appropriate input has to be
given in the calc: section, that is, layer composition and the respective calculators.
If energies & gradients are sent via sockets and the i-PI-protocol calc: can be
left empty but the appropriate socket addresses have to specified later in the opt:
section (see below). A separate socket/address is required for each layer.
Independent of the actual layer, pysisyphus always expects vectorial quantities with
appropriate shapes of the total system. The expected size of a force vector for a
system comprising N atoms is 3N. Vector components not belonging to frozen
atoms not belonging to a given layer are ignored and may be zero. It is expected,
that the full valid ONIOM forces for the total system are sent with the innermost layer.
These forces are then also used to check for optimization convergence.
A layered optimization is requested via type: layer in the opt: section. If
no further input is given, appropriate defaults will be used.
Micoiterations for the outer layers 0 to n-2 will be carried out using Cartesian
coordinates and regularized L-BFGS, while the innermost layer is optimized in internal
coordinates (type: tric) using rational function optimization (type: rfo).
If a more fine grained control is desired each layer can be controlled independently
by further keywords in the layers: list.
When layers: is present in the input all
layers must be specified, but empty list entries are possible (see the lines
only containing a comment in the example below).
One must start with the outmost layer; the last list item
corresponds to the innermost layer.
Geometry and/or optimizer setup of each layer is controlled by the geom: and opt:
sections in the list of layers with the usual keywords.
The usual keywords to control the optimization are supported, e.g., max_cycles,
thresh etc.
geom:
type:
cartesian
fn: [filename of initial geometry]
calc:
type:
oniom:
...
opt:
type:
layer:
[layers:
- # dummy for layer 0 (real system)
- # dummy for layer 1 layer
- # ...
# - # dummy for layer n-1 (model system)
- geom:
type:
cartesian:
opt:
type:
lbfgs:
]
# thresh: gau_loose
# max_cycles: 150
Below you can find a full example for the ONIOM2-optimization of Hexaphenylethane using pysisyphus' ONIOM implementation.
8.2. Example using pysisyphus' ONIOM
# 2-layer ONIOM optimization of Hexaphenylethan. The model
# system comprises the two carbons of the central ethan moiety.
# The model system is described GFN2-XTB while the real system is
# treated using GFN-FF.
geom:
type: cart
fn: |
68
C 1.2639 0.6380 0.4785
C 2.1638 1.9009 0.1349
C 2.4167 2.3403 -1.2006
C 3.2096 3.4576 -1.4755
C 2.7467 2.7036 1.1595
C 3.5368 3.8191 0.8739
C 3.7741 4.1922 -0.4416
H 2.0100 1.8572 -2.0700
H 3.3815 3.7580 -2.5014
H 2.5944 2.5089 2.1979
H 3.9618 4.4009 1.6820
H 4.3848 5.0587 -0.6598
C 0.7891 0.0052 -0.8980
C -1.0360 -0.5225 -2.5243
C -0.0432 -1.1375 -3.4197
C 1.2448 -1.1452 -3.0690
C -0.6577 0.0144 -1.3528
H -1.4492 0.4671 -0.7765
H -2.0793 -0.4955 -2.8238
H -0.3570 -1.5703 -4.3642
C 1.6732 -0.5526 -1.7810
H 1.9838 -1.5872 -3.7297
H 2.7287 -0.5813 -1.6050
C -0.0387 1.1726 1.2081
C -0.2929 2.5599 1.4341
C -1.0635 0.2863 1.6524
C -2.2174 0.7426 2.2933
C -2.4082 2.0998 2.5118
C -1.4505 3.0060 2.0772
H 0.3724 3.3423 1.1178
H -1.0153 -0.7691 1.5014
H -2.9711 0.0356 2.6161
H -1.6065 4.0659 2.2335
H -3.3046 2.4506 3.0066
C 2.1016 -0.4883 1.4312
C 3.4003 -1.0269 0.6933
C 4.4239 -0.1425 0.2420
C 5.5730 -0.6012 -0.4056
C 3.6509 -2.4147 0.4669
C 4.8038 -2.8631 -0.1829
C 5.7602 -1.9589 -0.6243
H 4.3778 0.9134 0.3928
H 6.3259 0.1044 -0.7337
H 2.9864 -3.1958 0.7882
H 4.9572 -3.9235 -0.3393
H 6.6530 -2.3115 -1.1243
C 2.5704 0.1468 2.8077
C 4.3159 0.7126 4.4630
C 3.3714 1.2844 5.3047
C 2.0334 1.2840 4.9354
C 3.9303 0.1495 3.2436
H 4.7333 -0.2864 2.6778
H 5.3575 0.6988 4.7582
H 3.6749 1.7164 6.2495
C 1.6376 0.7248 3.7185
H 1.2932 1.7141 5.5981
H 0.5868 0.7426 3.5290
C 1.2016 -1.7498 1.7785
C 0.9537 -2.1874 3.1154
C 0.6140 -2.5534 0.7571
C -0.1758 -3.6680 1.0474
C -0.4081 -4.0392 2.3643
C 0.1611 -3.3038 3.3950
H 1.3648 -1.7032 3.9822
H 0.7614 -2.3600 -0.2822
H -0.6047 -4.2504 0.2418
H -0.0069 -3.6027 4.4220
H -1.0186 -4.9050 2.5861
calc:
type: oniom
calcs:
real:
type: xtb
gfn: ff
pal: 6
high:
type: xtb
pal: 6
models:
# The "real" does not have to be described, only the "real" system.
high:
# Use ethan-carbons for the model system
inds: [0,34]
# Use appropriate calculator for model system.
calc: high
opt:
# Use LayerOpt
type: layers
thresh: gau
# When the 'layer' key is present, the user has to provide one
# entry per layer in the optimization. When the layer is empty default
# values are used. The layers are given from the outmost to the innermost
# layer.
# In the example above no other arguments for the outmost layer (# real dummy),
# while redundant internal coordinates are used for the innermost layer.
#
# When pysisyphus own ONIOM implementation is used no additional atom indices
# must be given here.
layers:
- # real dummy
- geom:
type: redund
do_hess: True
8.3. Example with sockets & i-PI-protocol
Whereas pysisyphus can figure out the layer composition for the microcycles when its own
ONIOM calculator is used, the user has to specify it when using sockets and the i-PI-protocol.
When covalent bonds between layers exist link atoms (LAs) must be introduced. Given a
layer with index n, atoms in layer n-1 that are replaced by LAs in the given layer
are called link atom hosts (LAH). A given layer and its LAHs in higher layers must
be optimized together.
Fig. 8.2 Layer structure of a ONIOM2 optimization using microiterations with link atoms. A given layer as well as connected link atom hosts in higher layers must be optimized simultaneously.
A simple example for a ONIOM2 optimization with microiterations is found below. Here, ethanal is optimized, with the model system comprising the carbonyl group (atoms 4, 5 and 6). This introduces a link atom between the two carbons 4 and 0, with carbon 0 becoming a LAH.
Fig. 8.3 ONIOM2 optimization of ethanal. All atoms surrounded by the dashed blue line are optimized in the innermost layer. Hydrogens are shown white, carbons grey and oxygen red.
In contrast to the first example using the native ONIOM implementation the user
can omit any input in the calc: section. Now the socket addresses have to given
for every layer, starting with the total system. For the total system the atom indices can
be omitted, as it is assumed that it comprises all atoms. For the remaining layers,
the indices of all moving atmos including LAHs in higher layers have to be given.
geom:
type: cartesian
fn: lib:acetaldehyd_oniom.xyz
calc:
opt:
type: layers
layers:
# Atom indices of the total system / (first layer) don't have to be
# specified.
- address: ./layer0sock
# Optimize atoms of the model system [4, 5, 6] and the link atom host 0.
- indices: [0,4,5,6]
address: ./layer1sock
thresh: gau
# The assert: section can be left out in actual inputs,
# it is here for testing purposes.
assert:
opt_geom.energy: -115.53522653
Another sensible choice for optimizing outer layers besides (regularized) L-BFGS may be preconditioned L-BFGS (type: plbfgs).
- class pysisyphus.optimizers.LayerOpt.LayerOpt(geometry, layers=None, **kwargs)[source]
Bases:
Optimizer- check_convergence(*args, **kwargs)[source]
Check if we must use the model optimizer to signal convergence.
- property layer_num: int
- property model_opt
Return the persistent optimizer belonging to the model system. We don't have to supply any coordinates to the optimizer of the most expensive layer, as it is persistent, as well as the associated geometry.
- class pysisyphus.optimizers.LayerOpt.Layers(geometry, opt_thresh, layers=None)[source]
Bases:
object
- class pysisyphus.optimizers.LBFGS.LBFGS(geometry, keep_last=7, beta=1, max_step=0.2, double_damp=True, gamma_mult=False, line_search=False, mu_reg=None, max_mu_reg_adaptions=10, control_step=True, **kwargs)[source]
Bases:
Optimizer- __init__(geometry, keep_last=7, beta=1, max_step=0.2, double_damp=True, gamma_mult=False, line_search=False, mu_reg=None, max_mu_reg_adaptions=10, control_step=True, **kwargs)[source]
Limited-memory BFGS optimizer.
See [1] Nocedal, Wright - Numerical Optimization, 2006 for a general discussion of LBFGS. See pysisyphus.optimizers.hessian_updates for the references related to double damping and pysisyphus.optimizers.closures for references related to regularized LBFGS.
- Parameters:
geometry (
Geometry) -- Geometry to be optimized.keep_last (
int, default:7) -- History size. Keep last 'keep_last' steps and gradient differences.beta (
float, default:1) -- Force constant β in -(H + βI)⁻¹g.max_step (
float, default:0.2) -- Upper limit for the absolute component of the step vector in whatever unit the optimization is carried out.double_damp (
bool, default:True) -- Use double damping procedure to modify steps s and gradient differences y to ensure sy > 0.gamma_mult (
bool, default:False) -- Estimate β from previous cycle. Eq. (7.20) in [1]. See 'beta' argument.line_search (
bool, default:False) -- Enable implicit linesearches.mu_reg (
Optional[float], default:None) -- Initial guess for regularization constant in regularized LBFGS.max_mu_reg_adaptions (
int, default:10) -- Maximum number of trial steps in regularized LBFGS.control_step (
bool, default:True) -- Wheter to scale down the proposed step its biggest absolute component is equal to or below 'max_step'**kwargs -- Keyword arguments passed to the Optimizer baseclass.
- class pysisyphus.optimizers.PreconLBFGS.PreconLBFGS(geometry, alpha_init=1.0, history=7, precon=True, precon_update=1, precon_getter_update=None, precon_kind='full', max_step_element=None, line_search='armijo', c_stab=None, **kwargs)[source]
Bases:
Optimizer- __init__(geometry, alpha_init=1.0, history=7, precon=True, precon_update=1, precon_getter_update=None, precon_kind='full', max_step_element=None, line_search='armijo', c_stab=None, **kwargs)[source]
Preconditioned limited-memory BFGS optimizer.
See pysisyphus.optimizers.precon for related references.
- Parameters:
geometry (
Geometry) -- Geometry to be optimized.alpha_init (
float, default:1.0) -- Initial scaling factor for the first trial step in the excplicit line search.history (
int, default:7) -- History size. Keep last 'history' steps and gradient differences.precon (
bool, default:True) -- Wheter to use preconditioning or not.precon_update (
int, default:1) -- Recalculate preconditioner P in every n-th cycle with the same topology.precon_getter_update (
Optional[int], default:None) -- Recalculate topology for preconditioner P in every n-th cycle. It is usually sufficient to only determine the topology once at the beginning.precon_kind (
Literal['full','full_fast','bonds','bonds_bends'], default:'full') -- What types of primitive internal coordinates to consider in the preconditioner.max_step_element (
Optional[float], default:None) -- Maximum component of the absolute step vector when no line search is carried out.line_search (
Literal['armijo','armijo_fg','strong_wolfe','hz',None,False], default:'armijo') -- Whether to use explicit line searches and if so, which kind of line search.c_stab (
Optional[float], default:None) -- Regularization constant c in (H + cI)⁻¹ in atomic units.**kwargs -- Keyword arguments passed to the Optimizer baseclass.