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 the occupation numbers will be in the range 0-1.

Note

In case of fractional occupation (i.e. non-integer occupations 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
    C    0.0005234578    0.0001727093    0.0002762926    0.0002311556
    C    0.0000008845    0.0007601116    0.0000014245    0.0007601017
    C    0.0005218907    0.0001737959    0.0002777874    0.0002325467
    O    0.0004150496    0.0002459280    0.0000738370    0.0004012048
    O    0.0004146230    0.0002469354    0.0000735124    0.0004003429
    H    0.0000668939    0.0001815939    0.0001866700    0.0001726131
    H    0.0001069923    0.0001412011    0.0001620722    0.0001367296
    H    0.0001069626    0.0001408593    0.0001624397    0.0001365849
    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.