*****************
Subtractive QM/MM
*****************

.. py:class:: NLayerSubtractive

Subtractive (ONIOM-like) QM/MM calculations can be specified with the ``NLayerSubtractive`` class as follows::

  my_frag = Fragment(coords=...)

  my_qm: chemsh._QM = ...
  qm_atom_indices: List[int] = ...
  qm_layer: chemsh.Layer = Layer(my_qm, qm_atom_indices)

  my_mm: chemsh._MM = ...
  mm_atom_indices: List[int] = ...
  mm_layer: chemsh.Layer = Layer(my_mm, mm_atom_indices)

  my_oniom = NLayerSubtractive(frag=...,
                 layers=(qm_layer, mm_layer), 
                 ...)

where ``frag`` is a :py:class:`<class Fragment`> object, ``layers`` is a tuple of objects (:class:`Layer`) that define the layers from the innermost to the outermost. In the example above, mm_atom_indices (i.e. the outermost atom indices) are optional. The above example is generalisable to any number of layers (>1), by supplying ``layers=`` with a tuple containing the desired layers. 

.. _`QM`: qm.html

.. _`MM interfaces`: mm.html

Once the :py:class:`NLayerSubtractive` driver theory has been set up, it can be used in tasks like
any other theory object, e.g. for a single-point calculation::

  sp = SP(theory=my_oniom, gradients=True)

  sp.run()

NLayerSubtractive options
=========================

.. py:attribute:: embedding

   Allowed values:

   * ``'mechanical'``: (default) Mechanical embedding

   * Electrostatic embedding will be implemented in a future release

.. py:attribute:: layers

   `Tuple` of `Layers` objects each defining non-overlapping regions of the molecular system and the corresponding computational methodologies


Layer
=====

.. py:class:: Layer

.. py:attribute:: method: _Interfaces

   The QM or MM method to be used in for the computation of this layer (will also be used to build subtractive terms for all preceding layers).

.. py:attribute:: atom_indices: Union[NDArray[np.int32], Iterable[int]] 

    A list of integers or a numpy array of integers containing all the atomic indices to be included in the current layer.
    (Optional in outermost layer where all the remaining atomic indices in the fragment are assumed to define the layer).

.. py:attribute:: gradients_eval: Union[Literal["analytic"], Literal["finitediff"]]

   Choice of method by which to determine the gradients in this layer:

   * ``'analytic'``: (default) Use analytic forces provided by `method`.

   * ``'finitediff'`` Use finite-difference (implemented by ChemShell) to approximate the forces.  

.. py:attribute:: link_atoms: Union[List[str], None]

   A list of strings containing the desired link atom types (Must be defined for MM layers, optional for QM layers: defaults to hydrogen). 

.. py:attribute:: link_atom_charges: Union[List[float], None]

   A list of containing the charges of the link atoms (must be defined in MM layers and any layer where the link atoms are manually defined).
   For QM layers using the default parameter for `link_atoms` link atom charges are 0.