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
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.