Dipole moment, Polarizability, Bond orders, Orbital Energies
************************************************************

This page in the AMS manual describes the calculation of the dipole
moment, the polarizability, and bond orders.

   Properties
      DipoleMoment Yes/No
      Polarizability Yes/No
      BondOrders Yes/No
   End

Note that because these properties are tied to a particular point on
the potential energy surface, they are found on the engine output
files. Note also that the properties are not always calculated in
every PES point that the AMS driver visits during a calculation. By
default they are only calculated in *special* PES points, where the
definition of special depends on the task AMS is performing: For a
geometry optimization properties would for example only be calculated
at the final, converged geometry. This behavior can often be modified
by keywords special to the particular running task.


Charges, Dipole Moment, Polarizability
======================================

   Properties
      Charges Yes/No
      DipoleMoment Yes/No
      DipoleGradients Yes/No
      Polarizability Yes/No
   End

"Properties"
   "Charges"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests the engine to calculate the atomic charges.

   "DipoleMoment"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests the engine to calculate the electric dipole moment
         of the molecule. This can only be requested for non-periodic
         systems.

   "DipoleGradients"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests the engine to calculate the nuclear gradients of the
         electric dipole moment of the molecule. This can only be
         requested for non-periodic systems.

   "Polarizability"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests the engine to calculate the polarizability tensor of
         the system.


Bond orders & Molecule detection
================================

Many engines can determine bond orders between atoms. For engines
based on force fields, these might just be the bond orders used
internally by the force field, while for quantum mechanical engines
the bond orders are usually determined by analyzing the results of the
quantum mechanical calculation, e.g. the electronic density. We refer
users to the manuals of the respective engine for details.

If bond orders are requested, but the engine does not implement a
method for determining bond orders, a bond guessing algorithm will be
used instead (see also BondOrders block). The bond guessing algorithm
(which is also used by the graphical user interface to guess bonds)
uses inter-atomic distances and atomic valences to guess the bond
orders.

You can request bond orders to be calculated in the "Properties" input
block:

   Properties
      BondOrders Yes/No
   End

"Properties"
   "BondOrders"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests the engine to calculate bond orders.
         For MM engines these might just be the defined bond orders
         that go into the force-field, while for QM engines, this
         might trigger a bond order analysis based on the electronic
         structure. For engines that do not have a bond order analysis
         method, a bond guessing algorithm will be used. See also the
         input options in the BondOrders block.

More options can be specified in the "BondOrders" block:

   BondOrders
      Method [Engine | Guess | EngineWithGuessFallback]
   End

"BondOrders"
   "Method"
      Type:
         Multiple Choice

      Default value:
         EngineWithGuessFallback

      Options:
         [Engine, Guess, EngineWithGuessFallback]

      Description:
         How to compute the bond orders when they are requested via
         the 'Properties%BondOrders' key.
         'Engine': let the engine compute the bond orders. The
         specific method used to compute the bond orders depends on
         the engine selected, and it may be configurable in the
         engine's input. Note: the calculation may stop if the engine
         cannot compute bond orders.                      'Guess': Use
         a bond guessing algorithm based on the system's geometry.
         This is the same algorithm that is used by the Graphical User
         Interface to guess bonds.
         'EngineWithGuessFallback': let the engine compute the bond
         orders (same as in 'Engine' option) but if the engine did not
         produce any bond orders, use the bond guessing algorithm as a
         fallback option.

Based on the bond orders, the AMS driver can analyze the atomic
connectivity graph in order to determine which sets of atoms together
constitute molecules. This allows for example to monitor the changes
in molecular composition during a reactive molecular dynamics
calculation. For molecular dynamics calculations this option is
enabled by default. For other tasks the molecular analysis can
explicitly be requested in the "Properties" block.

   Properties
      Molecules Yes/No
   End

"Properties"
   "Molecules"
      Type:
         Bool

      Default value:
         No

      Description:
         Requests an analysis of the molecular components of a system,
         based on the bond orders calculated by the engine.

Details of the molecule detection are configured in a dedicated block:

   Molecules
      AdsorptionSupportRegion string
      BondOrderCutoff float
   End

"Molecules"
   Type:
      Block

   Description:
      Configures details of the molecular composition analysis enabled
      by the Properties%Molecules block.

   "AdsorptionSupportRegion"
      Type:
         String

      GUI name:
         Adsorption support region

      Description:
         Select region that will represent a support for adsorption
         analysis. Adsorbed molecules will receive an '(ads)' suffix
         after name of the element bonded to the support. Such
         elements will be listed separate from atoms of the same
         element not bonded to the support, for example, HOH(ads) for
         a water molecule bonded to a surface via one of its H atoms.

   "BondOrderCutoff"
      Type:
         Float

      Default value:
         0.5

      Description:
         Bond order cutoff for analysis of the molecular composition.
         Bonds with bond order smaller than this value are neglected
         when determining the molecular composition.


Bond guessing algorithm
=======================

The algorithm to determine bond orders purely from atomic positions,
which is used if "BondOrders%Method" = "Guess", is roughly as follows:

* Each chemical element has a maximal number of bonds ("connectors"
  property taken from the PLAMS PeriodicTable class) and a covalent
  radius.

* Calculate interatomic distances

* Determine bond orders between all non-metal atoms

* Add bonds to any isolated H atoms (these bonds do not count towards
  the maximal number of bonds for an element).

* Find central atoms of common anions like carbonate, nitrate,
  sulfate, phosphate, and arsenate

* Add bonds between metal atoms and all other atoms, except the ones
  from the previous step

* Delete metal-metal bonds and metal-hydrogen bonds if the metal is
  bonded to enough electronegative atoms and not enough metal atoms.
  This means that the metal is likely a cation.

The calculated bonds are useful for visualization, but not necessarily
for setting up ForceField simulations.

The bond orders to metal atoms are always 1.

Note:

  The problem of finding molecular bonds for a given set of atoms in
  space does not have a general solution, especially considering the
  fact the chemical bond in itself is not a precisely defined concept.
  For every method, no matter how sophisticated, there will always be
  corner cases for which the method produces disputable results.
  Moreover, depending on the context (area of application) the desired
  solution for a particular geometry may vary. The algorithm used here
  gives good results for geometries that are not very far from the
  optimal geometry, especially consisting of lighter atoms. All kinds
  of organic molecules, including aromatic ones, usually work very
  well. Problematic results can emerge for transition metal complexes,
  transition states, incomplete molecules etc.

Warning:

  This method works reliably only for geometries representing complete
  molecules. If some atoms are missing (for example, a protein without
  hydrogens) the resulting set of bonds would usually contain more
  bonds or bonds with higher order than expected.

A version of this algorithm is also implemented in the PLAMS
Molecule.guess_bonds() method.


Orbital energies, HOMO-LUMO gap
===============================

Ab-initio or semiempirical engines, such as ADF or DFTB, internally
compute molecular orbital energies and occupations (see Summary of
engine capabilities). This information is printed to the AMSResults
section of the text output and is stored in the AMSResults section of
the binary engine results.

Example of molecular orbitals information printed to the text output:

   -----------------------
   Molecular orbitals info
   -----------------------
   n orbitals:                8
   n spin:                    1

   HOMO (eV):             -13.5919
   LUMO (eV):              -4.2036
   HOMO-LUMO gap (eV):      9.3883

      Index   Energy (eV)    Occupation
          8       12.3038         0.000
          7        9.4801         0.000
          6       -1.3196         0.000
          5       -4.2036         0.000 < LUMO
          4      -13.5919         2.000 < HOMO
          3      -14.8182         2.000
          2      -16.6543         2.000
          1      -20.6183         2.000

"n spin": 1 in case of spin-restricted or spin-orbit calculations, 2
in case of spin-unrestricted calculations.

"Occupations": in case of spin-unrestricted calculations the
occupation numbers will be in the range 0-2. In case of spin-
restricted calculation or spin-orbit calculations, the occupation
numbers will be in the range 0-1.

Note:

  In case of **fractional occupation** (i.e. non-integer occupation
  numbers), the definition of HOMO and LUMO is somewhat ill-defined.If
  fractional occupations are detected, a message will be printed to
  the text output.The logical variable
  "AMSResults%fractionalOccupation" in the engine .rkf results file
  indicates whether fractional occupation was detected.Here, we use
  the following convention:

  * The HOMO is the highest-energy orbital for which
    occupation[iOrbital] > occupied_threshold. The LUMO is HOMO+1.

  * if occupation[LUMO] < unoccupied_threshold we have integer
    occupations.

  * if occupation[LUMO] > unoccupied_threshold we have we have
    **fractional occupations**

  where occupied_threshold = 0.9 * max_occupation,
  unoccupied_threshold = 0.1 * max_occupation (max occupation = 2 for
  spin-restricted calculations and 1 for spin-unrestricted or spin-
  orbit calculations)

In case of spin-unrestricted calculation, the "smallest HOMO-LUMO gap"
is the smallest HOMO-LUMO gap irrespective of spin (i.e. min(LUMO) -
max(HOMO)).

The molecular orbital information is generally automatically computed
by engines supporting it.

You can explicitly request this in the input via the following
keyword:

   Properties
      OrbitalsInfo Yes/No
   End

"Properties"
   "OrbitalsInfo"
      Type:
         Bool

      Default value:
         No

      Description:
         Basic molecular orbitals information: orbital energies,
         occupations, HOMO, LUMO and HOMO-LUMO gap.

If you explicitly request "OrbitalsInfo" to be computed, but the
engine cannot compute it, AMS will stop with an error.

Molecular orbital information can only calculated in case of non-
periodic systems.


Uncertainties
=============

Some engines may report uncertainties for some properties. Currently
only the Hybrid engine supports uncertainties and only when running in
committee mode.

Uncertainties are reported to text output:

   Engine energy uncertainty (hartree)              0.01253685

   Uncertainty in Gradients (hartree/bohr)
     Index   Atom            d/dx            d/dy            d/dz       Magnitude
         1      C    0.0005234578    0.0001727093    0.0002762926    0.0002311556
         2      C    0.0000008845    0.0007601116    0.0000014245    0.0007601017
         3      C    0.0005218907    0.0001737959    0.0002777874    0.0002325467
         4      O    0.0004150496    0.0002459280    0.0000738370    0.0004012048
         5      O    0.0004146230    0.0002469354    0.0000735124    0.0004003429
         6      H    0.0000668939    0.0001815939    0.0001866700    0.0001726131
         7      H    0.0001069923    0.0001412011    0.0001620722    0.0001367296
         8      H    0.0001069626    0.0001408593    0.0001624397    0.0001365849
         9      H    0.0000668212    0.0001810509    0.0001866413    0.0001723906

Additionally, they can be obtained from the engine rkf file in the
"AMSResults%EnergyUncertainty", "AMSResults%GradientsUncertainty" for
the uncertainty in the energy in the energy and gradients
respectively. For convenience,
"AMSResults&GradientsMagnitudeUncertainty" is provided for
convenience, which reports the uncertainty in the magnitude of the
gradients based on the variance formula, or error propagation.

During "Task MolecularDynamics", the uncertainties are additionally
provided in the *History* section of *ams.rkf* in the variables
"EngineEnergyU(#)", "EngineGradientsU(#)" and
"EngineGradientsNormU(#)" (Magnitude).

It is possible to define an exit condition when uncertainties are too
high.

See also:

  In PLAMS, you can use dedicated getters to extract this information
  from AMSResults.
