Part 2: Spin-polarized Simulation
For the first part we carried out a structure optimization for the water molecule, H\(_2\)O. Here, we will run a structure optimization for the oxygen molecule, O\(_2\).
The contrast of the two examples illustrates a key lesson: The spin state of a molecule or solid matters. H\(_2\)O is a closed-shell molecule. However, O\(_2\) is actually spin-polarized in its gas-phase ground state. It has two unpaired electrons, which assume the same spin state to form a paramagnetic molecule in nature.
It is critical to consider the spin state of a given molecule or solid as an input parameter to an electronic structure code, on equal footing with the initial geometry guess. The spin state is NOT an automatic output of a simulation - already because spin-polarized systems usually have more than one spin state that could be stabilized, even if only one such spin state is the ground state and all other possible spin states are not. Click on the following information box for some further comments.
Non-spinpolarized vs Spin-polarized Systems
The H\(_2\)O molecule can safely be assumed to be non-spinpolarized in almost any practical scenario. Therefore, the simulation can be run in a nonmagnetic variant, without considering spin. For a molecule or solid that is truly non-magnetic, running a spin-polarized simulation would cost significantly more time and might even lead to an unstable self-consistent field cycle and/or a wrong spin state as the outcome.
In contrast, two possible spin states can be stabilized for O\(_2\): non-magnetic, i.e., spin-paired (higher energy and chemically much more reactive) or paramagnetic, i.e., two unpaired spins (lower energy and chemically much less reactive). In order to simulate the correct, physically desired state (usually, the ground state), it is essential to provide the relevant initial information to FHI-aims and, at the end of the simulation, verify that the desired state was in fact obtained.
Consider that you are likely sitting in front of a keyboard as you read this, immersed in an atmosphere with a significant partial pressure of the second-most chemically reactive element in the world: Oxygen. If this oxygen were non-magnetic, none of us would be sitting here in our present form. The fact that O\(_2\) is paramagnetic (spin-polarized) in its ground state is a critical determinant of its properties. No simulation will get the Earth-atmospheric, gas-phase properties of O\(_2\) right if spin is not considered.
The spin keyword
The spin-state of a material or molecule is controlled by keyword spin in FHI-aims.
For systems known to be non-spinpolarized, spin none is the correct keyword. This choice avoids unnecessary computational complexity and traps. We implicitly used this choice for H\(_2\)O, since spin none is the default if the spin keyword was not specified in control.in.
For suspected spin polarized systems, on the other hand, the true ground-state can only be obtained with a spin-polarized simulation and spin collinear is the correct keyword to be used. There is, however, an unavoidable complication for spin-polarized systems. Usually, more than one spin state can be stabilized in a simulation (and often even in experiment). Therefore, a spin-polarized simulation requires a correct initialization of the spin configuration of the material.
For spin-polarized simulations, a correct spin initialization must therefore also be provided as input to the simulation. This should be done in geometry.in, by adding specific keywords initial_moment after each atom that should carry an initial spin. For instance,
initial_moment 1.0
Why the initial_moment keyword(s) should be used in geometry.in
There are good reasons why we recommend against a blanket spin-polarized initialization of all atoms with a certain default initial moment. For any system of reasonable complexity, such an initialization will almost certainly correspond to a physically unreasonable initialization. For instance, a high-spin state (all spins oriented equally) will invariably be assumed, ruling out any antiferromagnetic solutions practically a priori.
Additionally, a blanket spin-polarized initialization can cause severe or insurmountable problems with convergence of the self-consistent field cycle, if the assumed spin configuration is far from the actual ground state. This can cost lots of computer time, spent on exploring entirely irrelevant parts of the electronic configurational space of a material or molecule. Thinking about the desired spin state a priori is a critical part of ensuring efficient simulations that rapidly lead to physically correct conclusions.
After a spin-polarized simulation is complete, the predicted spin state may be different from the state provided in the initialization. This final spin-state must be carefully checked at the end of a simulation in order to correctly verify the final result of a calculation. For many systems with stable spin ground-states, this task is straightforward. However, for systems with competing spin states, much more care and attention can be required to get the spin state right.
We therefore always recommend (in practice, we insist) on providing a thoughtful spin initialization in the geometry.in file of a spin-polarized simulation, using initial_moment for each atom that might carry a spin. There are usually multiple possible solutions and even severe technical pitfalls associated with incorrect spin guesses. An electronic structure code, in practice, cannot be able to guess on its own which one of potentially many possible physical scenarios a user is trying to explore.
Relaxation of the O\(_2\) Molecule
The relaxation process for O\(_2\) is very similar to the preceding example (H\(_2\)O), with the critical difference that the ground state of O\(_2\) is known to be spin-polarized. If we just followed the example of H\(_2\)O above, we would in fact obtain a perfectly well converged solution, but for the wrong spin state (zero spin-polarization) - an unphysical state of the molecule. Clearly, attention must be paid to how the calculation is initialized in order to arrive at the correct result.
The key difference is the use of spin collinear in control.in, and the specification of one or more initial spin moments for the atoms in geometry.in. The beginning of your control.in file should look as follows:
xc pbe
relativistic atomic_zora scalar
relax_geometry bfgs 5e-3
spin collinear
Also, we add a light basis set to control.in. This can be done by using the command
cat species_defaults/defaults_2020/light/08_O_default >> control.in
In order to pick a good initial spin moment, we do indeed require some chemical knowledge. Alternatively, we could do a systematic search by testing all reasonable combinations of positive and negative initial moments (remember, spin configurations can be parallel or antiparallel in a simple scalar-relativistic picture). With the systematic search, we would find all possible spin-polarized solutions of the O\(_2\) molecule. However, a look into the chemistry textbook convinces us that the O\(_2\) molecule in our atmosphere has two unpaired electrons, that is, one unpaired electron per O atom. This is why we choose to set initial_moment 1.0 in geometry.in of the initial structure. In this example, the initial distance between the O atoms is larger than the experimental value. We expect that the structure relaxation results in a more physical bond length of the O dimer.
atom 0.0 0.0 0.0 O
initial_moment 1.0
initial_charge 0.0
atom 0.0 0.0 1.9 O
initial_moment 1.0
initial_charge 0.0
You may note that we here also set use the keyword initial_charge to set the initial charge of each O atom to 0.0. Although this is not necessary for the O\(_2\) molecule, the example illustrates that an initial charge for each atom can also be specified in the geometry.in file. In certain cases, such as in some ionic materials, choosing the right initial charge can drastically improve the convergence of the self-consistent field cycle. For an example, please see the tutorial about simulating transition metal oxides with FHI-aims.
Once the input files are ready, begin the simulation as before by calling the FHI-aims executable (again, using the mpirun command) in your example folder.
Checking the results
After the relaxation is finished, the final structure of the molecule looks can be found in geometry.in.next_step and should look like this:
#
# This is the geometry file that corresponds to the current relaxation step.
# If you do not want this file to be written, set the "write_restart_geometry" flag to .false.
# aims_uuid : A9F3BD59-D533-47AE-A828-209E94D8C593
#
atom -0.00000000 -0.00000000 0.33687125 O
atom -0.00000000 -0.00000000 1.56312875 O
#
# What follows is the current estimated Hessian matrix constructed by the BFGS algorithm.
# This is NOT the true Hessian matrix of the system.
# If you do not want this information here, switch it off using the "hessian_to_restart_geometry" keyword.
#
trust_radius 0.2000000030
hessian_file
In the aims.out file resulting from this run, we can confirm that the FINAL total spin moment printed for the O\(_2\) molecule is indeed 2:
[...]
Current spin moment of the entire structure :
| N = N_up - N_down : 2.00
| S : 1.00
| J : 3.00
Highest occupied state (VBM) at -6.90564459 eV
| Occupation number: 1.00000000
| Spin channel: 1
Lowest unoccupied state (CBM) at -4.63949797 eV
| Occupation number: 0.00000000
| Spin channel: 2
Overall HOMO-LUMO gap: 2.26614662 eV.
| Chemical Potential : -5.44135698 eV
| Note that, for insulating systems, the printed 'chemical potential' value is not uniquely defined.
| It can be anywhere in the energy gap, as long as it correctly separates occupied and unoccupied states.
| In systems with a gap, the physically relevant chemical potential is the VBM or HOMO.
Self-consistency cycle converged.
[...]
aims.out to the Output Analyzer App of GIMS.
It is important to note that the "chemical potential" mentioned here refers to the approximate Fermi level between highest occupied and lowest unoccupied energy states. The O\(_2\) molecule has a non-zero HOMO-LUMO gap. From a computational point of view, the electronic "chemical potential" is therefore somewhat arbitrary since all that is needed is to place it somewhere in the HOMO-LUMO gap. For systems with a HOMO-LUMO gap (also semiconductors and insulators), the printed electronic "chemical potential" is therefore not unique and very meaningful, and it can fluctuate seemingly randomly within the gap, between different systems. Therefore, for systems with an energy gap, the highest occupied state (HOMO or VBM) would be a more reproducible choice for the reference energy of electrons.
So far, we have learned about how to obtain the local-minimum energy structure and total energy of simple molecules such as H\(_2\)O and O\(_2\). We also know about the importance of spin (charge) initialization for convergence of spin-polarized systems. We are now set to move on to the next example.
Solutions
You find all the solution to all the above exercises by clicking on the button below.