Electronic Configuration
************************

The next few keys can be used to specify the electronic configuration.
If you don't specify any such keys, certain defaults will apply. In
principle, the program will (by default) attempt to find the lowest-
energy spin-restricted (one-determinant) state. If SCF convergence is
problematic, the program may converge to an excited state, by which
(in this context) we mean a one-determinant state with a higher energy
than some other one-determinant state with the same net spin
polarization. In worse cases, the program may fail to converge to any
state at all. It is good practice to *always* verify which
configuration you actually have computed.

When you specify a particular configuration and/or net charge and/or
net spin-polarization of the system, the program will attempt to
compute that state accordingly, even if the requested configuration
has no physical or chemical meaning. The program has no knowledge
about the existence of materials and will simply try to carry out what
you tell it to do.


Charge and Spin
===============


Spin: restricted vs. unrestricted
---------------------------------

   Unrestricted Yes/No

"Unrestricted"
   Type:
      Bool

   Default value:
      No

   Description:
      By default, a spin-restricted calculation is performed where the
      spin alpha and spin beta orbitals are spatially the same.

   SpinPolarization float

"SpinPolarization"
   Type:
      Float

   Default value:
      0.0

   Description:
      The spin polarization of the system, which is the number of
      spin-alpha electrons in excess of spin-beta electrons.
      Specification is only meaningful in a spin-unrestricted
      calculation. However, specification is not meaningful in an
      unrestricted Spin-Orbit coupled calculation using the
      (non-)collinear approximation.

In an "unrestricted" calculation, spin-\alpha and spin-\beta MOs may
be spatially different and may have different occupation numbers. In
the case of spin-orbit coupling, "unrestricted" means that Kramers'
symmetry does not have to be satisfied. In that case, the key
"UNRESTRICTED" should be used in combination with the subkey
"SpinOrbitMagnetization NONCOLLINEAR" or "SpinOrbitMagnetization
COLLINEAR" of the key "Relativity".

The unrestricted mode roughly doubles the computational effort. The
actual numbers of spin-\alpha and spin-\beta electrons are controlled
by the keys "SpinPolarization" and OCCUPATIONS. Note carefully that
using *only* the keyword "unrestricted", without either
"SpinPolarization" or "Occupations" (or both), will not result in any
spin polarization. This implies that you would effectively perform a
spin-restricted calculation, but with increased computational effort.
Therefore, the program will check that in an unrestricted calculation
at least one of the keys "SpinPolarization" or "Occupations" is
applied.

The unrestricted feature is equivalent, in *ab initio* terminology, to
(spin-)unrestricted Hartree-Fock (UHF); the N-particle wave function
is a single determinant and not necessarily an eigenfunction of the
spin operator S^2.

A *restricted* calculation implies that the (spatial) orbitals *and*
the occupation numbers are identical for spin-\alpha and spin-\beta.

An unrestricted calculation with "SpinPolarization 0.0" (for instance
by not specifying the spin polarization at all) is, if one does not
break the spin symmetry, in fact a restricted run: it should give
exactly the same result as the restricted calculation, but it will use
more CPU time. If one does break the spin symmetry, for example with
the key "MODIFYSTARTPOTENTIAL" or the "SPINFLIP" option in the key
RESTART, the solution may also be a broken spin symmetry solution. For
example, you may want to start a calculation in broken symmetry with
spin-\alpha density on one fragment and spin-\beta density on another,
e.g. in a spin-unrestricted calculation of H_2 at large separation.

The Fock operator, both in an unrestricted and in a restricted run,
commutes with the spin operator Sz, but not (unless accidentally) with
S^2. The obtained one-determinant wave function may for instance be a
mixture of a singlet and a triplet state.

In an unrestricted calculation the expectation value of S^2 is
computed in ADF (note 29 in Ref. [1]). The evaluation of S^2 is not
entirely straightforward. DFT is essentially a one-particle formalism,
so the S-operator for the n-particle system has to be written out in
single-particle operators [2]. The equations used in ADF to calculate
the expectation value of S^2 can be found in Szabo and Ostlund [3].
Note that the so-called exact value (S_exact )^2, which is printed in
the ADF output, is defined as (S_exact )^2 = (|N_a - N_b |/2)(|N_a -
N_b |/2+1), where N_a and N_b are the number of spin-\alpha and
spin-\beta electrons, respectively. The expectation value of S^2 is
not calculated in a spin-orbit-coupled calculation.

If one includes the UNRESTRICTEDFRAGMENTS key, molecules that have
been calculated using the spin-unrestricted formalism can be employed
as fragments. In this case, the calculation on the full complex also
needs to use the spin-unrestricted formalism. See also the examples
Example: unrestricted fragments: CH3I and Example: unrestricted EDA:
Cu(C2H4)2.

A fair approximation to a computation with spin-unrestricted fragments
can be achieved with the key FRAGOCCUPATIONS, which uses spin-
restricted fragments, but allows for spin-unrestricted fragment
occupations. See also the examples Example: Spin-unrestricted
Fragments: H2, Example: Bond Energy analysis meta-GGA, (meta-)hybrids:
Zn2, Cr2, CrH and Example: Bond Energy analysis open-shell fragments:
PCCP.


Spin: unrestricted vs. restricted open-shell
--------------------------------------------

In unrestricted calculations the spin-\alpha and spin-\beta orbitals,
energies and occupations may differ. In restricted open shell
calculations the spin-\alpha and spin-\beta orbitals and energies
should be identical, only the occupations may differ. In ADF2023, a
restricted open-shell method was implemented for high-spin open-shell
molecules (ROSCF), similar to Ref. [11], in which case the one-
determinant wave function is an eigenfunction of S_z and S^2. The
implementation of the method is a work in progress. At the moment, it
is only valid for single-point calculations. This method cannot be
used with spin-orbit coupling or for many other properties. In
ADF2025, ROSCF can be used to calculate excitation energies of open-
shell molecules, including a method to include spin-orbit coupling
perturbatively.

The method implemented in ADF needs integer occupation numbers and a
positive spin polarization. It requires the subkey ROSCF of the key
SCF, the keyword UNRESTRICTED, and the key SPINPOLARIZATION. Part of
the input may look like:

   Occupations IntegerAufbau
   Unrestricted Yes
   SpinPolarization 2
   SCF
       ROSCF
       End
   End

The ROSCF options are:

   SCF
      ROSCF
         Alpha float_list
      End
   End

"SCF"
   "ROSCF"
      Type:
         Block

      Description:
         Settings for the ROSCF method.

      "Alpha"
         Type:
            Float List

         Default value:
            [0.5, 0.5, 0.5]

         Description:
            Coefficients to build the alpha-spin orbital contribution
            to the diagonal closed-, open-, and virtual-shell blocks
            of the Fock matrix. The beta-spin orbital contributions
            are 1.0 minus the alpha ones.


Unrestricted and Spin-Orbit Coupling
------------------------------------

In the case of spin-orbit coupling, there are two ways to perform
spin-polarized calculations, either using the collinear approximation
or the non-collinear approximation [4] [5]. Using the unrestricted
feature to assign different numbers of electrons to \alpha and \beta
spin cannot be applied as such when spin-orbit coupling is included,
since the electrons are not directly associated with spin-\alpha and
spin-\beta. For the collinear and non-collinear approximations, you
should use "Symmetry NOSYM" (see Symmetry key), and each level can
allocate 1 electron. Do not specify the spin polarization.

**Collinear approximation**

   Unrestricted Yes
   Symmetry NOSYM
   Relativity
      Level Spin-Orbit
      SpinOrbitMagnetization COLLINEAR
   End

The "SpinOrbitMagnetization" subkey is only relevant in the case of
Spin-Orbit coupling.

In the collinear approximation, at each point in space the spin
polarization has the same direction (default is in the direction of
the z-axis). Kramers' symmetry does not have to be satisfied. Symmetry
should be NOSYM. The default direction of the spin polarization can be
overridden using as argument "COLLINEARX" for spin polarization only
in the direction of the x-axis, or as argument "COLLINEARY" for spin
polarization only in the direction of the y-axis.

**Noncollinear approximation**

   Unrestricted Yes
   Symmetry NOSYM
   Relativity
      Level Spin-Orbit
      SpinOrbitMagnetization NONCOLLINEAR
   End

In the non-collinear approximation, at each point in space the spin
polarization can have a different direction. Kramers' symmetry does
not have to be satisfied. Symmetry should be NOSYM.

**Spin-orbit coupled spin-polarized start-up potential**

In the collinear and non-collinear approach, by default ADF uses an
artificial magnetic field in the z-direction in the first SCF cycle
that only affects spin (spin-Zeeman term). With the
"SpinOrbitMagnetization" key you can use an artificial magnetic field
that can have a different orientation in different regions. This is
similar to using a spin-polarized start-up potential, however, the
usage is different and with the key "SpinOrbitMagnetization" you can
also choose an arbitrary magnetization direction. Note that a
different orientation is typically only useful in the non-collinear
approach.

Remark: keys "MODIFYSTARTPOTENTIAL" and "RESTART%spinflip" cannot be
used in spin-orbit coupled calculations. The key
"SpinOrbitMagnetization" should be used instead.

   SpinOrbitMagnetization
      PerRegion
         Direction float_list
         Region string
         Strength float
      End
      Strength float
   End

"SpinOrbitMagnetization"
   Type:
      Block

   Description:
      Starting artificial spin-orbit magnetization at first SCF cycle

   "PerRegion"
      Type:
         Block

      Recurring:
         True

      Description:
         Defines the starting spin-polarization direction for all
         atoms in a region. Note that if this keyword is used multiple
         times, the chosen regions may not overlap.

      "Direction"
         Type:
            Float List

         Default value:
            [0.0, 0.0, 1.0]

         Description:
            The starting spin-polarization direction for a region.

      "Region"
         Type:
            String

         Description:
            The identifier of the region for which to define the
            starting spin-polarization direction.

      "Strength"
         Type:
            Float

         Default value:
            0.2

         Unit:
            Hartree

         Description:
            Strength of artificial spin-polarization in this region

   "Strength"
      Type:
         Float

      Default value:
         0.2

      Unit:
         Hartree

      Description:
         Strength of artificial spin-polarization at first SCF cycle


Net Charge
----------

The net charge of the molecule can be controlled with the key CHARGE
in the System block in the AMS part of the input:

   System header
      Charge float
   End

"System"
   "Charge"
      Type:
         Float

      Default value:
         0.0

      GUI name:
         Total charge

      Description:
         The system's total charge in atomic units.


Orbital occupations: electronic configuration, excited states
=============================================================

With the keys "OCCUPATIONS" and "IRREPOCCUPATIONS" you can specify in
detail the assignment of electrons to MOs.


Aufbau, smearing, freezing
--------------------------

   OCCUPATIONS Options

"Options"
   May contain IntegerAufbau, Keeporbitals, ElectronicTemperature,
   Smear, Freeze, or Steep:

   "IntegerAufbau"
      Electrons are assigned to MOs according to the Aufbau principle,
      trying to use integer occupations, also in case of degeneracy at
      the Fermi level. ADF normally may use fractional occupation
      numbers in that case. Note that for multi-dimensional irreps,
      using the subkey IntegerAufbau may still not prevent fractional
      occupation numbers from being used.

   "Keeporbitals=NKeep"
      Until SCF cycle Nkeep electrons are assigned to MOs according to
      the Aufbau principle, using at each cycle the then current
      orbital energies of the MOs. Thereafter the KeepOrbitals feature
      is applied. As soon as this is activated the program will on
      successive SCF cycles assign electrons to the MOs that maximally
      resemble - in spatial form - those that were occupied in a
      'reference cycle number'. The default for Nkeep is 20, except:

      1. When orbital occupations for MOs are specified explicitly in
         the data block of the occupations key, these apply
         throughout.

      2. In a Create run fixed occupations are derived from a database
         in the program.

      3. When electron smearing is explicitly turned on by the user
         (see the Smear option below) Nkeep is by default 1,000,000 so
         the program will 'never' compare the spatial forms of MOs to
         determine the occupation numbers. The 'reference cycle
         number' is by default the previous cycle, which will suppress
         jumps in the spatial occupations during the SCF development
         while, on the other hand, allowing the system to let the
         more-or-less-frozen configuration relax to self-consistency.

   "Freeze"
      Occurrence of this word in the option list specifies that the
      'reference cycle number' will be the cycle number on which the
      KeepOrbitals feature is activated: during all subsequent SCF
      cycles the program will assign electrons to MOs that resemble
      the MOs of that specific SCF cycle. This may be used when the
      MOs of that cycle are already reasonably close to the final
      ones, and it will suppress unwanted step-by-step charge-
      transfers from occupied to empty orbitals that are very close in
      energy. By default this option is not active.

   "ElectronicTemperature=T"
      Invokes orbital population using the Fermi-Dirac distribution
      with temperature "T" (in Kelvin). This may be used to achieve
      SCF convergence in an otherwise problematically converging
      system. Note that the default SCF acceleration method A-DIIS is
      not compatible with non-aufbau orbital occupations so it is
      disabled automatically and the classical Pulay DIIS scheme is
      used. For this reason it is strongly recommended to use either
      MESA or one of the methods from the LIST family in combination
      with "ElectronicTemperature" or "Smear".

   "Smear=Smear1[,Smear2,Smear3,...,Smear10]"
      Another form of specifying the electronic temperature (in
      Hartree). The smearing parameter Smear*N* is related to the
      "ElectronicTemperature" in such a way that the occupation number
      for any orbital outside the smearing interval will deviate from
      an integer value (0 and 1 for spin-unrestricted, 0 and 2 for
      spin-restricted) by no more than 1% (0.01 electron for spin-
      unrestricted, 0.02 for spin-restricted). More precisely, a
      smearing value of 1 Hartree is equivalent to an electronic
      temperature of 34285 Kelvin. When a comma-delimited list of
      values is specified, the next value from the list is picked
      after SCF has converged, and the SCF is continued. This way one
      can specify a list of gradually decreasing values to get a sort
      of annealing effect. NOTE: No spaces are allowed when specifying
      the list of values for Smear.

   "Steep=Lambda[,Nmax]"
      The occupation number for each orbital is updated according to a
      steepest-descent method (Ref: F. W. Averill and G. S. Painter,
      Phys. Rev. B **46**, 2498 (1992)). During an SCF cycle, the
      occupation number for each new orbital is initially determined
      by decomposing the old charge density with new orbitals. Then,
      the occupation numbers are modified so that the total energy of
      the system decreases. The Lambda parameter gives the coefficient
      for the charge transfer in 1/au units. The second parameter,
      Nmax, is an additional limit for the amount of charge transfer.
      Nmax would be useful for early cycles when the Lambda parameter
      gives too large a charge transfer. Too small an Nmax results in
      irregular behavior in SCF convergence. In the case of difficult
      SCF convergence, you should make mixing and Lambda smaller. From
      our experience, Nmax=0.1 or 0.2 is usually OK. This method
      should be used with the DIIS method turned off (DIIS N=0), and
      the choice of the mixing parameter in the SCF cycle is also
      important. This option is especially useful for systems with
      many quasi-degenerate orbitals around the Fermi level. For
      instance, cluster models of surface systems usually suffer from
      dangling bonds and should be converged with this method. Note
      though that slow convergence is an intrinsic feature of this
      method, so one should specify a large limit for the number of
      SCF cycles, say 500 or even 1000, depending on the cluster size.

   "OptimizeSpin=Delta", "OptimizeSpinRound=Delta"
      The spin polarization can be optimized by minimizing the
      difference between Fermi level values for spin-alpha and spin-
      beta MOs. This method works reliably only when a reasonably
      large finite electron temperature is used (300K is found to be
      sufficiently large for small transition metal clusters). The
      polarization value is optimized up to a specified tolerance
      Delta (fractional values are allowed) by step-wise reduction of
      the root-bracketing interval for the objective function
      f=FermiLevel(alpha)-FermiLevel(beta). The SCF procedure is
      converged for each spin polarization. When the interval becomes
      smaller than Delta the optimization is considered converged.
      After that the SCF may optionally be repeated with the spin
      polarization set to the nearest even or odd number (depending on
      the total number of electrons in the system) or, in other words,
      with the total number of electrons per spin rounded to the
      nearest integer. The method is invoked using either the
      OptimizeSpin or the OptimizeSpinRound Occupation option. As
      mentioned above, it must be used together with the
      ElectronicTemperature option. The OptimizeSpin and
      OptimizeSpinRound argument specifies the tolerance value. The
      OptimizeSpinRound option specifies that after polarization has
      converged the SCF must be completed with the total number of
      electrons per spin rounded to the nearest integer. It should be
      noted that, as usual, the Unrestricted and SpinPolarization keys
      must be specified as for any other spin-unrestricted
      calculation. The argument of the SpinPolarization key sets the
      initial spin polarization value from which the optimization is
      started. The result of the optimization may depend on the
      initial value because the objective function may have more than
      one root.

Notes about the occupations options:

* When occupation numbers are explicitly defined via the block
  "IRREPOCCUPATIONS" (see next section), the Smear option cannot be
  used.

* The Aufbau principle does not determine or adjust the distribution
  of electrons over spin-\alpha versus spin-\beta in an unrestricted
  calculation. This aspect is controlled by the key "SpinPolarization"
  and by any explicit occupations in the data block of occupations.

* Smearing cannot be used in combination with the KeepOrbitals option.


Explicit occupation numbers
---------------------------

   IrrepOccupations
     irrep orbitalnumbers
     irrep orbitalnumbers
     ...
   End

"irrep"
   The name of one of the irreducible representations (not a
   subspecies) of the point group of the system. See the Symmetry for
   the irrep names as they are used in ADF.

"orbitalnumbers"
   A series of one or more numbers: the occupation numbers for one-
   electron *valence* orbitals in that irrep. The orbitals are ordered
   according to their energy eigenvalue; higher states than those
   listed get an occupation number of zero.

   For degenerate representations such as the 2-dimensional
   E-representations or the 3-dimensional T-representations, you must
   give the *total* occupation, i.e. the sum over the partner
   representations; ADF assigns each partner an occupation equal to
   the appropriate fraction of what appears here.

   In an unrestricted calculation, two sequences of numbers must be
   specified for each irrep; the sequences are separated by a double
   slash (//). The first set of numbers is assigned to the spin-\alpha
   orbitals, the second set to the spin-\beta orbitals. Example
   unrestricted calculation in symmetry NOSYM with two unpaired
   electrons:

      IrrepOccupations
         A 28 // 26
      End
      SpinPolarization 2
      Symmetry NOSYM

   Note that this is not meaningful in an unrestricted Spin-Orbit
   coupled calculation using the (non-)collinear approximation, where
   one should use one sequence of occupation numbers for each irrep.

Notes about the occupations data block:

* When specifying electron configurations, all valence electrons in
  the calculation must be explicitly assigned to MOs and the
  IRREPOCCUPATIONS keyword must be used. In this context the concept
  *valence electrons* and hence *valence orbitals* is not necessarily
  identical to what you may normally assume to be the valence space of
  an atom or molecule. The meaning of *valence* is here strictly
  defined as whatever electrons are outside the frozen core. It
  depends therefore on the level of frozen core approximation applied
  in the calculation. This traces back to the Create runs in which the
  basic atoms were generated that are now used to build the molecule.

* When for some irrep there is a rather long list of occupation
  numbers, corresponding to *consecutive fully occupied* states, you
  can combine these numbers and enter their sum instead: ADF knows the
  maximum occupation for an irrep, and when you put a larger number
  the program will split it up. For instance, if you give for the
  *p*-representation (in a single atom calculation):

     P 17 3

  ADF will interpret this as:

     P 6 6 5 3

  i.e. the occupation number 17 is interpreted as denoting two fully
  occupied p-shells and the remaining five electrons in the next
  higher shell. This example also illustrates how to specify an
  excited state: here we have defined a hole in the third p-shell.

* Fractional occupation numbers in input are allowed. For a discussion
  of the interpretation of fractional occupation numbers see Ref. [6].
  The program even allows you (technically) to use a non-integer total
  number of electrons, whatever the physical meaning of such a
  calculation is.

* The data block of occupations is not parsed (see the section Input).
  The program does not replace expressions by their value and it does
  not recognize constants or functions defined with the define key.

* In a numerical frequencies run (without symmetric displacements) the
  symmetry used internally in the program is NOSYM, irrespective of
  any Schönfliess symbol in the input file. As a consequence the
  program will recognize only the A representation (the only irrep in
  nosym), but not the representations belonging to the input point
  group symmetry. (The symmetry in the equilibrium geometry, defined
  by the input Schönfliess symbol, is used to enhance efficiency and
  stability in the construction of the matrix of Force constants).


CHARGE and SPINPOLARIZATION vs. IRREPOCCUPATIONS
------------------------------------------------

The contents of the data block of IRREPOCCUPATIONS, if used, defines
the total number of valence electrons and hence the net total charge.
In an unrestricted run it also defines the net spin polarization. If
the keys "CHARGE" and "SPINPOLARIZATION" are also used, the program
will check that the specifications are consistent.

We strongly recommend that you always specify the net total charge and
spin polarization with "CHARGE" whenever explicit occupation numbers
are supplied with "IRREPOCCUPATIONS", so that the program can check
that your occupation numbers result in the total charge and spin
polarization you intend.


Create mode
-----------

In Create mode occupation numbers are predefined (see Appendix
Elements of the Periodic Table), and these are applied unless you
specify occupations in the input yourself. Conceivably this may result
in a non-Aufbau configuration. In Create mode the program always
operates as if the occupations were set in the input.


Frozen core approximation
=========================

**Frozen core vs. pseudopotentials**

Pseudopotentials are not supported. The frozen core approximation is
automatic in a normal (Fragment mode) calculation and is defined by
the basic atomic fragments. The data file used in the Create run
specifies the frozen core for the atom, which is then used in all
molecules that incorporate that atomic fragment.

**Core Potentials**

In the standard approach the Coulomb potential and the charge density
due to the atomic frozen core are computed from the frozen one-
electron orbitals. ADF stores the computed core density and core
potential for each atom type in the molecule on a file TAPE12.
Alternatively, you may attach a file with (core) potentials and
densities. The file must have the same structure as the standard
TAPE12. It should contain one or more sections, each with the core
information for one type of atom. With the block COREPOTENTIALS you
specify the core file and (optionally) which sections pertain to the
distinct atom types in the molecule.

   COREPOTENTIALS corefile
     {atomtype index}
     {atomtype index}
     ...
   end

"corefile"
   The file with core potentials and charge densities. The name may
   contain a path.

"atomtype"
   One of the atom type names as defined by atoms.

"index"
   Points to the core section on the attached file that applies to the
   atom type. Different atom types may use the same section. A non-
   positive index tells the program that the atoms of that type don't
   have a frozen core. If the information on the corresponding
   fragment file (or data file in Create mode) indicates the contrary
   the program will abort with an error message.

If the key is used as a simple key (specifying only the core file) the
sections on the file are associated with the atom types in order: the
first section is used for the first atom type, et cetera. This is
overruled by applying the block form. However, since the key *must*
have the core file as argument, the block form requires that you apply
the continuation symbol: an ampersand ("&"), separated from the core
file name by at least one blank.

If you omit an atom type from the data block it gets a zero index (no
core).

The attached file may contain more sections than used in the
calculation, and the indices specified in the data block don't have to
be in ascending order, consecutive, or cover a specific interval.

When a file with non-standard (e.g. relativistic) cores is attached
and used in the calculation of an atom or molecule, and the result is
used as a fragment in a subsequent calculation, you should attach and
use the same core potentials again. Otherwise, the program will
internally compute the standard core potentials and hence implicitly
employ a different fragment than you may think, i.e. a fragment with
other properties. ADF will not check anything in this respect and
COREPOTENTIALS should therefore be handled with great care.

The primary application of the corepotentials option is to include
(scalar) relativistic corrections in the (frozen core part of the)
Fock operator. The relativistic core potentials can be computed with
the auxiliary program dirac (see the RELATIVITY keyword).


Spin-polarized start-up potential
=================================

The Fock matrix for the first SCF cycle is built by combining the Fock
matrices from the fragment files. When performing a Restart, the
initial Fock Matrix may be read from the specified "adf.rkf" file.


Spin-flip method for broken symmetries
--------------------------------------

It is possible to exchange alpha and beta electrons for selected atoms
when performing a restart from a previous spin-unrestricted
calculation. Note that for spin-orbit coupled calculations one should
look at the section on how to make a spin-orbit coupled spin-polarized
start-up potential.

In many cases, one wishes to perform a calculation of a low-spin
complex where spin-density is positive on some atoms and negative on
the others. It is usually very difficult to achieve SCF convergence if
one starts from scratch. Sometimes, the "ModifyStartPotential" feature
(see next section) helps with this problem but sometimes it does not.
A more robust way is to first perform a high-spin calculation and then
modify the resulting adf.rkf file by "flipping" the spin on some
atoms. This file then can be used to restart a subsequent low-spin
calculation.

Such a "flipping" can be performed during restart by specifying the
"SpinFlip" keyword in the Restart input block as shown below:

   Restart
      File HighSpin.results/adf.rkf
      ! SpinFlip keyword is followed by the indices of the
      ! atoms for which the flipping will be performed
      SpinFlip 1
   End

See also:

  The example Broken spin-symmetry: Fe4S4 and the AMS-GUI tutorial
  Spin Coupling in Fe4S4 Cluster demonstrate the use of the spin-flip
  option.


Modify the starting potential
-----------------------------

In some applications you may want to modify the initial potential
(from the restart file or the fragment files), see also the previous
section. This is achieved with the block "ModifyStartPotential". It
allows you to scale the potential so as to represent user-chosen
amounts of spin-\alpha and spin-\beta density on some or all of the
fragments. This will adjust the spin-\alpha and spin-\beta initial
potentials. In spin-orbit coupled calculations the
"ModifyStartPotential" key cannot be used; see instead the section on
how to make a spin-orbit coupled spin-polarized start-up potential.

The "ModifyStartPotential" option applies only to *unrestricted*
calculations, of course. It may be used to help the program find a
particular state. This might, for instance, be hard to find otherwise
due to the alpha-beta symmetry in the start-up situation. It may also
be useful to speed up the SCF convergence in case you know what the
final distribution of spin-\alpha and spin-\beta density over the
molecule will approximately be.

   ModifyStartPotential {specification}
   { frag alfa // beta
    frag alfa // beta
    ...}
   end

"specification"
   Must be *two numbers*, ASPIN and BSPIN, if provided at all. They
   specify the (relative) amounts of spin-\alpha and spin-\beta fit
   density to define the spin-dependent potential at the first SCF
   cycle. The coefficients retrieved from the fragment files (or from
   the restart file in case of a SCF restart) are scaled accordingly.
   This will not affect the *total* amount of fit density: the
   absolute values of ASPIN and BSPIN play no role, only their ratio.
   In case of a restart run the restart file must have been generated
   in a *restricted* calculation, while the continuation run must be
   an *unrestricted* one.

If no argument is given, a data block must be supplied with records
frag alfa // beta. This is very much similar to the main option with
ASPIN and BSPIN: you specify ASPIN and BSPIN now separately for each
fragment. This involves somewhat more input but increases the
possibilities for tuning the initial potential. Again, this can be
applied only in an unrestricted calculation. It cannot be used in a
restart: the affected fit coefficients are those from the fragment
files, while in an SCF restart run these are ignored and replaced by
the coefficients on the adf.rkf (TAPE21) restart file.

Each line specifies a frag with its corresponding ASPIN and BSPIN fit
partitioning. If frag is the name of a fragment *type*, the specified
ASPIN-BSPIN is applied to all individual fragments of that type.
Alternatively an *individual* fragment can be specified, using the
format fragtype/n, where *n* is an index between one and the total
number of fragments of that type. In such a case the ASPIN-BSPIN data
applies only to that particular fragment while different values may be
supplied for the other fragments of the same type.

It is allowed to specify for certain fragment types individual
fragments and for other fragment types only the type. Duplicate
specifications are not allowed; an individual fragment must not be
specified if its fragment type is also specified as a whole. If the
data block form is used, only the fit coefficients of the referenced
fragments are affected. For the unreferenced fragments the fit
densities are used as they are defined on the corresponding fragment
files.

The SCF convergence of a spin-unrestricted calculation usually
improves when you start with potentials that correspond to the correct
ratio of spin-\alpha and spin-\beta electrons. By default
ASPIN=BSPIN=0.5, as implied by the spin-restricted start density of
the fragments or restricted molecule.

The total amount of fit density used on the first iteration is defined
by the sum-of-fragment densities (or the density on the restart file).
This may be different from the total number of electrons in the actual
calculation. On the second SCF cycle the fit density will internally
be normalized so as to represent the correct number of electrons.

The block-form of the key makes the start-up of broken symmetry
calculations easy. For example, one may want to start a calculation in
broken symmetry with spin-\alpha density on one fragment and
spin-\beta density on another, e.g. in a spin-unrestricted calculation
of H_2  at large separation. It is particularly useful for larger
systems, e.g. for magnetic coupling between spin-polarized magnetic
centers, as in Fe-S complexes [7]: start with oppositely polarized Fe
centers, but with, for instance, the remaining bridge and terminal
ligands unpolarized.

See also:

  The example Localized Hole: N2+ demonstrates the use of the
  "ModifyStartPotential" option.


Unrestricted fragments
======================

In ADF you can use spin-unrestricted fragments if the key
"UNRESTRICTEDFRAGMENTS" is used.

   UNRESTRICTEDFRAGMENTS
   UNRESTRICTED

If spin-unrestricted fragments are used, the full complex must also be
calculated spin-unrestricted, which means that the key "UNRESTRICTED"
is required, and in addition the "SPINPOLARIZATION" key and/or the
"IRREPOCCUPATIONS" block key must be included. The key
"FRAGOCCUPATIONS" should not be included. The spin-\alpha and
spin-\beta occupations are read from the spin-unrestricted fragment
files. For spin-restricted fragments the spin-\alpha and spin-\beta
occupations will be set equal and sum to the spin-restricted
occupation numbers.

The unrestricted fragments can be self-consistent: different numbers
of spin-\alpha and spin-\beta electrons usually result in different
spatial orbitals and different energy eigenvalues for spin-\alpha and
spin-\beta when you go to self-consistency.


Prepared for bonding
--------------------

Typically an unrestricted electron configuration for the fragments is
used, such that the Pauli repulsion between the fragments is minimal,
which means that the fragments are so-called 'prepared for bonding'.
For example, if one has two fragments that both have one unpaired
electron, one would put the unpaired electron of the first fragment in
the spin-\alpha orbital and the unpaired electron of the second
fragment in the spin-\beta orbital. If one wants to calculate the
electron pair bonding separately, see key REMOVEALLFRAGVIRTUALS.


Simulated unrestricted fragments with key FRAGOCCUPATIONS
=========================================================

In the method that simulates unrestricted fragments with the key
FRAGOCCUPATIONS, the fragments from which the molecule is built must
be spin-restricted, that is, the fragment files must be result files
of spin-restricted calculations. For purposes of analysis, however, it
may be desirable in some applications to build your molecule from
fragments that have an *unrestricted* electronic configuration. This
can be simulated as follows.

You tell ADF that you want to *treat* the fragments as if they were
unrestricted; this causes the program to duplicate the one-electron
orbitals of the fragment: one set for spin-\alpha and one set for
spin-\beta. You can then specify occupation numbers for these spin-
unrestricted fragments, and occupy spin-\alpha orbitals differently
from spin-\beta orbitals.

Of course, the unrestricted fragments that you use in this way are not
self-consistent: different numbers of spin-\alpha and spin-\beta
electrons usually result in different spatial orbitals and different
energy eigenvalues for spin-\alpha and spin-\beta when you go to self-
consistency, while here you have spatially identical fragment
orbitals. Nevertheless it is often a fair approximation which gives
you a considerable extension of analysis possibilities.


Prepared for bonding
--------------------

Typically an unrestricted electron configuration for the fragments is
used, such that the Pauli repulsion between the fragments is minimal,
which means that the fragments are 'prepared for bonding'. For
example, if one has two fragments that both have one unpaired
electron, one would put the unpaired electron of the first fragment in
the spin-\alpha orbital and the unpaired electron of the second
fragment in the spin-\beta orbital. If one wants to calculate the
electron pair bonding separately, see key REMOVEALLFRAGVIRTUALS.


FRAGOCCUPATIONS
---------------

   FRAGOCCUPATIONS
    fragtype
      irrep spin-a // spin-b
      irrep spin-a // spin-b
      ...
    subend
    fragtype
      irrep spin-a // spin-b
      ...
    subend
   end

"fragtype"
   One of the fragment types and functions as a (block type) subkey.
   The data block for the subkey ends with the standard end code for
   block type subkeys (subend).

"irrep"
   One of the irreducible representations (irreps) for the point group
   symmetry that was used in the computation of that fragment.

"spin-a // spin-b"
   Two sequences of occupation numbers, which will be applied to the
   spin-\alpha and spin-\beta versions of the Fragment Orbitals. The
   sequences must be separated by a double slash (//).

The sum of spin-\alpha and spin-\beta occupations must, for each
fragment orbital in each irrep separately, be equal to the total
(restricted) occupation of that orbital as it is stored on the
fragment file. In other words: you can only change the distribution
over spin-\alpha and spin-\beta electrons within one orbital.

(Without this restriction the spatial distribution of the total (sum
over spins) fragment charge density would be changed, leading to an
incorrect bonding energy analysis after the calculation).

Be aware that in more-dimensional irreps (E, T) the number of
electrons in a fully occupied orbital is input as the dimension of the
irrep times the one-electron orbital occupation. Compare the key
IRREPOCCUPATIONS.

For irreps that are not mentioned in this input block, and hence for
all irreps of fragment types that are not mentioned at all, the
spin-\alpha and spin-\beta occupations will be set equal, which is, of
course, what they are on the (restricted) fragment file.

For an example of applying this option see Ref. [8].


Remove Fragment Orbitals
========================

By default all fragment orbitals (the MOs of the fragment
computation), which are stored on the fragment file, are used as basis
functions for the overall molecule. You can remove one or more of
these fragment orbitals from the basis set of the molecule. This may
be useful for special analyses, like calculating the electron pair
bonding in case one has open-shell fragments, or for instance to study
the effect of deleting all virtual MOs of a particular fragment (CSOV
analysis, Constrained Space Orbital Variation). It may also enhance
the efficiency since you effectively reduce the size of the basis set,
but you should be aware of the potential effects on the results.

The pure orbital interaction effect of forming electron pair bonding
between open-shell molecules can approximately be calculated with a
bond energy analysis in which all virtual orbitals are removed from
the fragments, see Ref. [9]. For calculating the effect of electron
pair bonding, it is best to specify an unrestricted electron
configuration for the fragments, such that the Pauli repulsion is
minimal. Removing all virtuals from an open-shell fragment that is
calculated spin-restricted means that all fragment orbitals with zero
occupation are removed. Thus, for example, a singly occupied fragment
orbital will not be removed. This singly occupied fragment orbital
will result in a spin-\alpha and a spin-\beta fragment orbital. In
combination with other singly occupied fragment orbitals they may form
an electron pair bonding combination, but also an antibonding
combination. In practice this means that the orbital interaction
calculated with a bond energy analysis in which all virtual orbitals
are removed from open-shell fragments might be due to more than
electron pair bonding. The situation is not much different if one uses
unrestricted fragments, because technically always the same number of
spin-\alpha and spin-\beta orbitals must be removed, and the removed
orbitals must all have zero occupation.

If one wants to remove all virtual fragment orbitals, use the key
REMOVEALLFRAGVIRTUALS.

   REMOVEALLFRAGVIRTUALS

If one wants to remove all virtual fragment orbitals except on the
first fragment, use the key REMOVEOTHERFRAGVIRTUALS.

   REMOVEOTHERFRAGVIRTUALS

If one does not want to remove all virtual fragment orbitals, then one
should use the block key REMOVEFRAGORBITALS.

   REMOVEFRAGORBITALS
    fragtype
      subspecies nremove
      subspecies nremove
      ...
    subend
    fragtype
      subspecies nremove
      ...
    subend
    ...
    (etc.)
    ...
   end

"fragtype"
   One of the fragment types in the system. Any subset of the
   available fragment types can be used here as subkey. The subkeys
   are block type keys; their data blocks end with subend.

"subspecies"
   One of the subspecies of the irreducible representations of the
   point group symmetry that was used in the calculation of the
   fragment itself. This requires of course that one knows the
   symmetry that has been used for the fragment calculation.

"nremove"
   The number of fragment orbitals of the pertaining representation
   that will not be used as basis functions for the overall system.
   The *highest* (in energy eigenvalue) nremove orbitals are
   discarded. You must not remove *occupied* fragment orbitals.

By default (omission of the key) all fragment orbitals are used in the
basis set for the system.

**Important Note**

It is imperative that any removal of fragment orbitals will not break
the symmetry of the molecule. This consideration is relevant when, for
instance, two different subspecies of a fragment irrep contribute to
different partner subspecies in one of the irreps of the molecule. In
such a case, when one removes an orbital in such a fragment
subspecies, its partner orbital should also be removed. If this is
violated an error may occur or the results will simply be wrong. Quite
likely, the program will detect the error, but this may occur only in
the final (analysis) stage of the calculation so that a lot of CPU
time may have been wasted.

Example: consider a single-atom fragment, computed in atom symmetry,
used as fragment in a c(lin) molecule and assume that the p:x and p:y
fragment orbitals contribute to respectively the pi:x and pi:y
subspecies of the molecule. Then, when you remove one or more p:x
fragment orbitals, you should also remove the same number of p:y
fragment orbitals. Practical cases may be more complicated and
whenever you use this key, make sure that you have fully analyzed and
understood how the fragment irreps combine into the molecular symmetry
representations. Hint: run the molecule, without removing any fragment
orbitals, and stop at an early stage after the program has computed
and printed the build-up of the molecular SFOs from the fragment
orbitals. To control early aborts via input, use the key STOPAFTER.


CDFT: Constrained Density Functional Theory
===========================================

CDFT is a tool for carrying out DFT calculations in the presence of a
constraint. The nature of the constraint is general in theory;
however, in the current implementation the user can constrain the
charges or the spins of a set of moieties (as identified by sets of
atoms) to be a specific real number given in the input. Note that the
use of CDFT as implemented in ADF is an **expert option**, and it is a
**work in progress**. Implemented in ADF by M. Pavanello and P. Ramos
[11], based on the method described in Ref. [10]. At the moment
"SYMMETRY NOSYM" and an all-electron basis set are required.

Warning:

  CDFT is an expert option that may not work for your system or the
  specific constraint you had in mind.

The simplest way to run CDFT is using the following keyword
combination

   CDFT
     NCONSTRAINTS 1
     NATOMSPERSET N
     THEATOMS atom1 atom2 ... atomN
     CONSTRAINTS charge
   END

All the CDFT block key options are:

   CDFT
      AllAtoms Yes/No
      AnalyticalHessian integer
      ChargeAndSpin Yes/No
      Constraints float_list
      DoNotOptimize Yes/No
      ExcitedCDFT Yes/No
      InitialMultipliers float_list
      MaxIter integer
      Metric Yes/No
      NAtomsPerSet integer_list
      NConstraints integer
      OnlyCharge Yes/No
      OnlySpin Yes/No
      PopType [yukawalike | fuzzyvoronoibecke | fuzzyvoronoifermi]
      Print [low | medium | high]
      SelfConsistent Yes/No
      StepSize float
      TheAtoms integer_list
      Threshold float
   End

"CDFT"
   Type:
      Block

   Description:
      CDFT is a tool for carrying out DFT calculations in the presence
      of a constraint.

   "AllAtoms"
      Type:
         Bool

      Default value:
         No

      Description:
         If AllAtoms is true, then TheAtoms is overridden and all the
         atoms in the active fragment are included in the set.

   "AnalyticalHessian"
      Type:
         Integer

      Default value:
         0

      Description:
         This will calculate the analytical derivative of the energy
         w.r.t. the Lagrange multiplier up to the specified SCF
         iteration. This key is not recommended due to the high
         computational cost that comes with it. The calculation is
         equivalent to a ground state Hessian, and it is carried out
         with the full sum-over-states formula.

   "ChargeAndSpin"
      Type:
         Bool

      Default value:
         No

      Description:
         will constrain both the charge and the spin

   "Constraints"
      Type:
         Float List

      Description:
         The values of the constraints. If CHARGEANDSPIN, constraints
         to the alpha and beta electrons need to be specified
         sequentially. One more electron => CONSTRAINTS -1.0. One less
         electron => CONSTRAINTS 1.0. If the CDFT type is EXCITEDCDFT,
         CONSTRAINTS=1.0 is recommended. Other values are technically
         possible but have not been tested yet.

   "DoNotOptimize"
      Type:
         Bool

      Default value:
         No

      Description:
         If true, the multipliers chosen in INITIALMULTIPLIERS will
         not be optimized and will be constant throughout the entire
         SCF procedure.

   "ExcitedCDFT"
      Type:
         Bool

      Default value:
         No

      Description:
         will generate an excited state with CONSTRAINTS number of
         ALPHA electrons constrained to occupy the virtual space of a
         ground state reference calculation. This is the essence of
         the eXcited Constrained DFT (XCDFT) method(P. Ramos, M.
         Pavanello, Low-lying excited states by constrained DFT,
         Journal of Chemical Physics 148, 144103 (2018)
         https://doi.org/10.1063/1.5018615) for the calculation of
         low-lying single excitations. XCDFT is found to correctly
         reproduce the energy surface topology at conical
         intersections between the ground state and the first singly
         excited state and can also accounts for the condensed phase
         effects in solvated chromophores where typical Delta SCF
         methods variationally collapse.

   "InitialMultipliers"
      Type:
         Float List

      Description:
         If available, a guess for the Lagrange multipliers can be
         entered.

   "MaxIter"
      Type:
         Integer

      Default value:
         200

      Description:
         Maximum number of CDFT iterations. CDFT carries out a loop
         nested inside the SCF cycle.

   "Metric"
      Type:
         Bool

      Default value:
         No

      Description:
         Relevant for XCDFT. In the XCDFT method orthogonality is not
         imposed between the KS-orbitals of ground and excited states.
         If METRIC is specified, the degree of mixing of the single
         excited state with the ground state or high-order excitations
         is calculated. Three parameters are calculated: p, m and d.
         The parameters p and m will give information about the amount
         of mixing with the ground state, while parameter d will
         determine the mixing with high order excitations. Additional
         information about the origin of these parameters can be found
         in the literature (P. Ramos, M. Pavanello, Low-lying excited
         states by constrained DFT, Journal of Chemical Physics 148,
         144103 (2018) https://doi.org/10.1063/1.5018615)

   "NAtomsPerSet"
      Type:
         Integer List

      Description:
         The number of atoms in each moiety (set).

   "NConstraints"
      Type:
         Integer

      Default value:
         1

      Description:
         This specifies the number of sets of atoms to be considered.
         For example, if the user wishes to constrain a positive
         charge on one part of the system, and a negative charge on
         another part, NCONSTRAINTS should be set to two. There is no
         limit on the number of constraints. However, SCF convergence
         becomes an issue with more than 2 constraints. Note:
         NCONSTRAINTS>1 is untested.

   "OnlyCharge"
      Type:
         Bool

      Default value:
         Yes

      Description:
         Will constrain only the charge, letting spin relax (and
         potentially delocalize)

   "OnlySpin"
      Type:
         Bool

      Default value:
         No

      Description:
         Will constrain only the spin

   "PopType"
      Type:
         Multiple Choice

      Default value:
         yukawalike

      Options:
         [yukawalike, fuzzyvoronoibecke, fuzzyvoronoifermi]

      Description:
         The population analysis chosen for determining the
         constraint.

   "Print"
      Type:
         Multiple Choice

      Default value:
         low

      Options:
         [low, medium, high]

      Description:
         Print level and debugging.

   "SelfConsistent"
      Type:
         Bool

      Default value:
         No

      Description:
         Self-Consistent CDFT

   "StepSize"
      Type:
         Float

      Default value:
         0.5

      Description:
         The amount of the Lagrange multipliers step taken in each
         CDFT iteration

   "TheAtoms"
      Type:
         Integer List

      Description:
         The atom numbers of the moieties in the input geometry order.
         If NCONSTRAINTS is larger than 1, the sets of atoms are
         entered as a single list.

   "Threshold"
      Type:
         Float

      Default value:
         1e-10

      Description:
         The threshold for convergence of the CDFT constraints. The
         tighter the SCF convergence criteria, the tighter the
         THRESHOLD should be.


FOCDFT: Fragment Orbital Constrained DFT
========================================

ADF has a fragment-oriented approach. FOCDFT is a fragment-orbital-
based method for constraining (integer) charges and/or spins on
fragments. In practice this method works best if there are no covalent
bonds, or if one does not want covalent bonds between the fragments.
FOCDFT uses the so-called Löwdin orbitals constructed by ADF at the
start of an SCF calculation on the whole complex. These Löwdin
orbitals are reasonably localized, such that one can assign a fragment
to each Löwdin orbital. Next, these Löwdin orbitals are used to
constrain charges and/or spins on fragments.

Symmetry can be used in FOCDFT, although one typically should not use
symmetry-equivalent fragments. FOCDFT gradients cannot be calculated.
FOCDFT is not implemented for spin-orbit coupled calculations.

To better understand the method, it is relevant to see how the Löwdin
orbitals are constructed in ADF. The fragment MOs are called fragment
orbitals (FO) in ADF. If there are symmetry-equivalent fragments ADF
generates symmetry combinations of the FOs and uses the symmetrized
fragment orbitals (SFOs) as basis functions. In case of FOCDFT one
typically should not use symmetry-equivalent fragments, which means
that in that case there is no distinction between SFOs and FOs. The
so-called Löwdin orbitals are Löwdin-orthonormalized core-
orthogonalized fragment orbitals. This Löwdin orthonormalization is
done in several steps. In case of frozen core calculations the
fragment orbitals are first orthogonalized on the frozen core orbitals
on the other fragment(s). In case of an all electron calculation this
step is skipped. Next the (core orthogonalized) occupied fragment
orbitals are orthonormalized. The occupied-occupied part of the
overlap matrix in terms of (core orthogonalized) occupied fragment
orbitals is diagonalized. The results are used to do a Löwdin
orthonormalization (S^{-1/2}) of the (core orthogonalized) occupied
fragment orbitals. Next the (core orthogonalized) valence fragment
orbitals are Gram-Schmidt orthogonalized on the Löwdin-orthonormalized
(core orthogonalized) occupied fragment orbitals, and the overlap
between these functions is calculated and these virtual orbitals are
then Löwdin-orthonormalized. In practice the Löwdin orbitals coming
from the Löwdin-orthonormalized (core orthogonalized) occupied
fragment orbitals resemble the occupied fragment orbitals quite well.
The other Löwdin orbitals (that have more virtual fragment character)
may differ more from the virtual fragment orbitals they originate
from, especially if the basis set is larger. Note that the number of
virtual fragment orbitals increases if the basis set increases,
whereas the number of occupied fragment orbitals remains constant.

In FOCDFT the Fock matrix is calculated on the basis of the Löwdin
orbitals. The off-diagonal part of the Fock matrix between Löwdin
orbitals on one fragment and Löwdin orbitals on another fragment is
set to zero, such that a block diagonal Fock matrix is constructed.
Due to this block structure each eigenfunction of the Fock matrix will
be a linear combination of Löwdin orbitals on one fragment. This means
that the eigenfunctions are automatically localized on a fragment. By
occupying the orbitals with the lowest eigenvalues on each fragment,
it is very easy to constrain (integer) charge and/or spin on such a
fragment.

Usage:

   FOCDFT
      Charge float_list
      SpinPolarization float_list
   End

"FOCDFT"
   Type:
      Block

   Description:
      Localize orbitals on fragments and constrain charge and spin
      polarization on fragments.

   "Charge"
      Type:
         Float List

      Description:
         Constrain charge on fragments. Number of charges should be
         the same as number of fragments, and sum should add up to
         total charge.

   "SpinPolarization"
      Type:
         Float List

      Description:
         Constrain spin polarization on fragments. Number of spin
         polarization values should be the same as number of
         fragments, and sum should add up to total spin polarization.
         Calculation should be spin unrestricted.

The list of fragment charges and spin should be given in the same
order as the fragments in the FRAGMENTS key. Especially if one
includes charge constraints, the resulting electron configuration may
not be Aufbau.

FOCDFT can be used in the calculation of charge transfer integrals and
exciton transfer integrals.

[1] R.E. Bulo, A.W. Ehlers, S. Grimme and K. Lammertsma,
    *Vinylphosphirane.Phospholene Rearrangements: Pericyclic
    [1,3]-Sigmatropic Shifts or Not?* Journal of the American Chemical
    Society 124, 13903 (2002)

[2] R. Pauncz, *Spin Eigenfunctions*, ISBN13: 9780306401411, 1979, New
    York: Plenum Press

[3] A. Szabo and N.S. Ostlund, *Modern Quantum Chemistry*, ISBN13:
    9780070627390, 1st ed. revised ed. 1989: McGraw-Hill

[4] H. Eschrig and V.D.P. Servedio, *Relativistic density functional
    approach to open shells*, Journal of Computational Chemistry 20,
    23 (1999)

[5] C. van Wüllen, *Spin densities in two-component relativistic
    density functional calculations: Noncollinear versus collinear
    approach*, Journal of Computational Chemistry 23, 779 (2002)

[6] S.G. Wang and W.H.E. Schwarz, *Simulation of nondynamical
    correlation in density functional calculations by the optimized
    fractional orbital occupation approach: Application to the
    potential energy surfaces of* O_3 and SO_2, Journal of Chemical
    Physics 105, 4641 (1996)

[7] L. Noodleman, and E.J. Baerends, *Electronic Structure, Magnetic
    Properties, ESR, and Optical Spectra for 2-Fe Ferredoxin Models by
    LCAO-Xa Valence Bond Theory*, Journal of the American Chemical
    Society 106, 2316 (1984)

[8] F.M. Bickelhaupt, N.M. Nibbering, E.M. van Wezenbeek and E.J.
    Baerends, *The Central Bond in the Three CN* Dimers NC_CN, CN-CN,
    and CN-NC: Electron Pair Bonding and Pauli Repulsion Effects*,
    Journal of Physical Chemistry 96, 4864 (1992)

[9] F.M. Bickelhaupt, M. Solà, C. Fonseca Guerra, *Highly polar bonds
    and the meaning of covalency and ionicity -- structure and bonding
    of alkali metal hydride oligomers*, Faraday Discussions 135, 451
    (2007)

[10] Q. Wu, T. Van Voorhis, *Direct optimization method to study
     constrained systems within density-functional theory*, Physical
     Review A 72, 024502 (2005)

[11] P. Ramos, M. Pavanello, *Constrained subsystem density functional
     theory*, Physical Chemistry Chemical Physics 18, 21172 (2016)

[12] P. Ramos, M. Pavanello, *Low-lying excited states by constrained
     DFT*, Journal of Chemical Physics 148, 144103 (2018)

[13] Z. Rinkevicius, I. Tunell, P. Sałek, O. Vahtras, and H. Ågren,
     *Restricted density functional theory of linear time-dependent
     properties in open-shell molecules*, Journal of Chemical Physics
     119, 34 (2003)
