:sequential_nav: next .. _tutorial-basic-basis-optimization: Basis sets - Optimization ========================= .. dropdown:: **Do you want to get all tutorial materials on your local machine?** :color: primary You can get all files and set up your local environment. :ref:`Do that now ` before proceeding. Using the default basis sets generated by SIESTA might be enough for some applications, but doing some degree of manual optimization of the basis sets may help to achieve better results with similar computational costs. This is especially advisable when dealing with extremely large systems, where going for high basis set cardinality (triple-zeta, quadruple-zeta) is not really an option. *Note that the concepts presented in this tutorial apply to both molecules and crystals. Don't hesitate to use them for your own systems, whatever they are!* .. dropdown:: Additional information :color: info As a general background for this topic, you can watch the 2021 talks by Emilio Artacho: * Part I: https://www.youtube.com/watch?v=puweAzW9ZNs * Part II: https://www.youtube.com/watch?v=AlY2SDb1ta0 .. button-ref:: tutorial-basis-set-optimization-start :color: success :outline: :expand: :shadow: Take me to the exercises! :bdg-info:`background` General Concepts --------------------------------------- Cardinality: Multiple-zeta and polarization ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ SIESTA uses atomic-centered basis sets, which are the product of a radial function (R(r)) and a spherical harmonic (Y(m,l)). The basis set *cardinality* is what determines the total amount of basis functions that we have per atom. In a minimal basis set, one would have one basis function per valence orbital. So for a Carbon atom, we would get one 2s function and three 2p functions. This is called a *single-zeta* (SZ) basis set, which contains one function per orbital. We can, of course, increase the quality of this basis set by increasing the amount of radial functions that describe the same orbital. So in a *double-zeta* basis set (DZ), we have two functions per orbital. This means that for C, we would now have eight functions instead of the minimal four. This can go up indefinitely, to triple-zeta (TZ), quadruple-zeta (QZ), and so on. The last part that contributes to cardinality are the *polarization orbitals*. These increase the "flexibility" of a basis set, providing and extra set of orbitals of a higher *l*. So in the case of a double-zeta *polarized* basis set (DZP), our C atom will now add 5 d-type functions to the mix (one *l* higher than the p-orbitals), thus increasing the total to 13 functions. Save for a few exceptions, a DZP basis is more than enough to get a good description. In general, **increasing the cardinality of a basis set greatly increases the cost**, since the most time-consuming part of SIESTA has a cost proportional to the *cube* of the amount of basis functions. Cut-off radii ^^^^^^^^^^^^^ Remember that a basis function in SIESTA are the product of a radial function and a spherical harmonic. The spherical harmonic is determined by the quantum numbers of the orbital; as such, it is not something we can really modify. However, we can do work on the radial part. **The cut-off radius of a basis function indicates the point after which the radial function becomes zero**. For multiple-zeta basis sets, the matching radius indicates the point after which the n-zeta orbitals coincide with the first-zeta. In general, the larger the cut-off radii, the better described will be our system. However, this also comes at a cost; in particular, an increase in computational costs. So how do we choose a cut-off radius that is large enough to ensure good quality results, but also does not increase our simulation times unnecessarily? Basis Enthalpy ^^^^^^^^^^^^^^ The basis enthalpy is a fictitious magnitude that adds an extra term to the energy, which depends on the volume of the orbitals; essentially, **this can tell us when the basis orbitals are getting needlessly large**. This extra term is essentially P*V, where V is the orbital volume and P is a fictitious "pressure". This pressure is an input we will use in SIESTA. **For bulk materials** such as metals (Cu, Fe) or semiconductors (Si, Ge), the default value of **0.2 GPa for the basis pressure produces reasonable results.** However, **for small atoms** in isolated molecules (C, O, H, etc.) this pressure results in extremely short cut-off radii. **One typically reduces the basis pressure to a tenth of the default value, i.e. 0.02 GPa.** **How do you know whether your orbitals are large enough?** A good rule of thumb is to consider that orbitals with a cut-off radius of less than 5.0 Bohr are very small, while a cut-off radius larger than 10.0 Bohr is unnecessarily large. Of course, there are a few exceptions: 1s and 2s orbitals of the first two rows of the periodic table (specially H) can be slightly smaller than 5.0 Bohr, but going below 4.0 Bohr should be avoided. In the oposite case, when dealing with surfaces or 2D materials one might need larger, diffuse orbitals in order to better describe the area around the exposed atoms: these diffuse orbitals might easily reach 10 to 12 Bohr. More on this topic is covered in :ref:`tutorial-basic-basis-extra`. Optimization and transferability ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ When optimizing basis sets, the concept of **transferability** is one of the most important ideas to keep in mind. You might get the ultimate basis set that gives the perfect results for a given set of atoms in a given set of coordinates; however, if even a slight re-arrangement of those atoms produces terrible results, that means your basis set is not transferable. Ideally, one should aim for highly transferable basis sets that also produce reasonable results. This ties in with a key aspect of basis set optimization: defining your optimization function, i.e., the magnitude that you will be using to monitor how well your basis is performing. One good magnitude to look at is the total energy of a given system, but even better is to use the concept of **basis enthalpy**. While it might be tempting to optimize basis sets using bond lengths, interaction energies, or even band structures, these typically result in poorly transferable basis sets. The use of these magnitudes as optimization functions is not encouraged; however, these can be used *after* the optimization is done in order to benchmark the behaviour of the newly optimized basis set. Having defined the optimization function, one needs to define which are the optimization variables, i.e. the things that we will modify manually to improve our basis sets. While SIESTA has a ton of possibilities to modify your basis set, the most important targets for optimization are the orbital radii: the cut-off radii for first-zeta orbitals, and the matching radii for multiple-zeta. Polarization orbitals can also be optimized separately, but the way these are generated in SIESTA is generally good and they usually do not require further optimization. .. note:: To summarize: * We aim for good-quality, highly-transferable basis sets. * Typically we will be modifying the cut-off radii of first-zeta and the matching radii of multiple-zeta orbitals. * We will take the basis enthalpy (or the total energy) of a system as our main way of optimizing the basis sets. * We need to have a way to test if our basis set is transferable enough. .. _tutorial-basis-set-optimization-start: :bdg-success:`basic` Rough optimization of the basis set -------------------------------------------------------- .. note:: This tutorial shows just **one possible method** to optimize your basis sets. Different people tend to use different techniques, **there is no "one true way"**. This method does, however, work reliably across different kinds of systems. .. note:: We will stick to double-zeta polarized (DZP) basis sets for this tutorial, which should be enough for most calculations with SIESTA. .. figure:: images/water.png :width: 300px :align: center :alt: (For SIESTA, the water molecule is three points in space.) For SIESTA, the water molecule is three points in space: one representing oxygen (red) and two representing hydrogen (white). .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *01-Eshift_SplitNorm* Molecular dynamics of water (or aqueous systems in general) is still one of the most common types of simulations. In this tutorial, we will optimize the basis set for a single water molecule, and we will use the energy of a water dimer as a way to benchmark the quality of our basis functions. You can find the input files for this first section here: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`water.fdf` .. dropdown:: A peek at the fdf file .. literalinclude:: work-files/01-Eshift_SplitNorm/water.fdf :linenos: It is also important to note that we have set up a relatively big simulation box. Since water is an isolated molecule and not a bulk metal, we want to avoid orbital overlaps between one molecule and its periodic images: .. literalinclude:: work-files/01-Eshift_SplitNorm/water.fdf :linenos: :lineno-match: :start-at: LatticeConstant :end-at: %endblock LatticeVectors Take a look at the following lines in the *water.fdf* file: .. literalinclude:: work-files/01-Eshift_SplitNorm/water.fdf :linenos: :lineno-match: :start-at: PAO.EnergyShift :end-at: PAO.SplitNorm These two options allow us a rough control of the quality of the basis set. :ref:`PAO.EnergyShift` allows us to control the cut-off radii of the first-zeta orbitals: the lower the energy shift, the larger the cut-off radii. :ref:`PAO.SplitNorm` allows us to control the matching radii of the multiple-zeta orbitals: the lower the value, the larger the matching radii (i.e., more similar to the first-zeta orbitals). The values presented in the input file are the defaults. We want orbitals that are good quality but not needlessly large, so we will be relying on the Basis Enthalpy (as defined above) to guide our optimization. Since we are using small atoms, the variable :ref:`BasisPressure` will be set to 0.02 GPa, as you can see in the following line: .. literalinclude:: work-files/01-Eshift_SplitNorm/water.fdf :linenos: :lineno-match: :start-at: BasisPressure :end-at: BasisPressure You can run this fdf as usual: .. code-block:: console $ siesta < water.fdf > water.out Now that we have all set, lets start by running a single-point calculation for different values of :ref:`PAO.EnergyShift`. We will compare the results for the basis enthalpy in each case. You can get the basis enthalpy from the general output under "(Free)E+ p_basis*V_orbitals", or from the first line of the *.BASIS_ENTHALPY* file: .. code-block:: console $ grep "E + p_basis" water.out .. dropdown:: Output .. code-block:: (Free)E + p_basis*V_orbitals = -482.064557 eV .. code-block:: console $ cat h2o.BASIS_ENTHALPY .. dropdown:: Output .. code-block:: Basis enthalpy [eV] : -482.06455744756806 Basis Harris enthalpy [eV] : -482.06455749949799 Free energy [eV] : -482.13276932085228 Harris energy [eV] : -482.13276937278221 Orbital volume [Ang**3] : 546.43734768656418 Average basis pressure [eV/Ang**3] : 1.2483018148921527E-004 Average basis pressure [GPa] : 2.0000000000000000E-002 *Test the following values for the energy shift 0.01 Ry, 0.001 Ry, and 0.0001 Ry. Compare the basis enthalpy for all three values. Which of the three performs better?*. .. dropdown:: Answer: Basis enthalpy for different values of energy shift :color: success The basis enthalpy for the three values of the energy shift are as follows:: 0.01 Ry -482.064557 0.001 Ry -482.435979 0.0001 Ry -482.344446 We can see that the basis enthalpy is minimized for the value of 0.001 Ry. 0.0001 Ry might provide a lower *total energy*, but the orbitals become too large. Now we will test different values for :ref:`PAO.SplitNorm`. *Set the energy shift to the one we found above, 0.001 Ry. Test the values 0.45, 0.55 and 0.65 for the Split Norm. Which of the three gives the better basis enthalpy?* .. dropdown:: Answer: Basis enthalpy for different values of split norm :color: success The basis enthalpy for the three values of the split norm are as follows:: 0.45 Ry -482.696472 0.55 Ry -482.741379 0.65 Ry -482.621525 The best value for the split norm in this case is 0.55 Ry. If you have a look at the *.BASIS_ENTHALPY* file, you will see that :ref:`PAO.SplitNorm` does not have a strong impact on the total orbital volume when compared to energy shift: what you are seeing here, essentially, is how the value of the matching radii affects the total energy of the system. .. note:: In real-life applications, you would typically optimize both the energy shift and the split norm iteratively, as they might not be totally independent of each other. Be sure to explore a larger array of values for both parameters, although differences below 1 meV in the basis enthalpy are usually not that relevant. .. dropdown:: Additional information: The BasisPressure.Specs block :color: info Starting from SIESTA 5.4, you can now specify a value for the basis pressure in the :ref:`BasisPressure.specs` block, which allows you to set different values for different atomic species. You can do it by using the atomic number or the species label. For example:: %block BasisPressure.Specs H 0.02 GPa # Hydrogen Z 8 0.0003 GPa # Oxygen %endblock BasisPressure.Specs This is useful for systems that have a mix of small and large atoms. .. dropdown:: Additional information: How do the orbital shapes look like? :color: info If you can plot the radial part of the basis orbitals generated by SIESTA, you can see the actual effect of :ref:`PAO.EnergyShift` and :ref:`PAO.SplitNorm`. Let's take for example the 2s orbital of oxygen. The first zeta is only affected by the energy shift, and, unless the change is very drastic, the overall shape of the orbital doesn't change much: .. image:: images/o2p-rough.png :width: 500px :align: center The second zeta, however, is affected by both the energy shift and the split norm, and the changes in shape are much more noticeable. Take into account that the actual physical orbital in SIESTA is taken as the difference from the first zeta, so the plot is not really :math:`\phi_2s(zeta=2)`, but ctually :math:`\phi_2s(zeta=2)-\phi_2s(zeta=1)`: .. image:: images/o2p2z-rough.png :width: 500px :align: center If you check the output files, you will see that the cut-off radius for the first zeta and the matching radius for the second zeta coincide with the point at which the radial function becomes zero. .. dropdown:: Additional information: older versions of SIESTA :color: info In SIESTA 5.4.0, the output for the basis enthalpy was changed to include additional information. For previous versions, the information on the general output has to be grepped like this:: $ grep "E+ p_basis" water.out In addition the *.BASIS_ENTHALPY* file used to be split in two different files: the BASIS_ENTHALPY file and the BASIS_HARRIS_ENTHALPY file. The contain the same information, but formatted differently: .. code-block:: -482.06455744612362 The above number is the electronic (free)energy: -482.13276931940788 Plus the pressure : 1.3595723686972974E-006 ( 2.0000000000000000E-002 GPa) times the orbital volume (in Bohr**3): 3687.5420025311782 Starting from SIESTA 5, we greatly improved the defaults for generating basis sets. If, for any reason, you want to reproduce these results with an older version, try adding :ref:`PAO.SoftDefault` and ``PAO.OldStylePolOrbs`` to your input file:: PAO.SoftDefault T PAO.OldStylePolOrbs F These will guarantee closer results, but there are more things involved. Check the manual and the corresponding release notes for SIESTA 5. .. :ref:`PAO.OldStylePolOrbs` is not included in .. the manual... :bdg-success:`basic` Calculating the binding energy of a water dimer -------------------------------------------------------------------- .. figure:: images/dimer.png :width: 300px :align: center :alt: (A water dimer) A water dimer, i.e. two water molecules sharing a hydrogen bond. .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *02-WaterDimer* Now is the time to put our optimized basis set to test. We will use as an example the binding energy of a water dimer, i.e.: :math:`E_{binding} = E_{dimer} - 2*E_{monomer}` The input files are the same as before, with *monomer.fdf* being the single water molecule we used before: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`dimer.fdf` and :download:`monomer.fdf` .. dropdown:: A peek at dimer.fdf .. literalinclude:: work-files/02-WaterDimer/dimer.fdf :linenos: You will notice that both files have the default values for the energy shift and split norm: .. literalinclude:: work-files/02-WaterDimer/monomer.fdf :linenos: :lineno-match: :start-at: PAO.EnergyShift :end-at: PAO.SplitNorm Run both files as usual: .. code-block:: console $ siesta < monomer.fdf > monomer.out $ siesta < dimer.fdf > dimer.out Take note of the total energy for both systems (remember to look for Total in the output file). Then, change the :ref:`PAO.SplitNorm` and :ref:`PAO.EnergyShift` values in the *monomer.fdf* and *dimer.fdf* to the values obtained in the previous section (if you didn't do it, pick 0.001 Ry for the energy shift and 0.55 for the split norm). Run the calculations again. *How different is our binding energy for the dimer?* .. dropdown:: Answer: Basis enthalpy for different values of split norm :color: success For the default basis set options, the total energy of the monomer is -482.136281 eV, while the energy of the dimer is -964.717861 eV. This gives a binding energy of -445 meV. For the optimized basis set, the total energy of the monomer is -482.884135 eV, and the energy of the dimer is -966.016972 eV. This gives a binding energy of -248 meV. You can see now that the binding energy with the lightly optimized basis setis significantly different. Does this mean that the default basis sets are awful? Not really, because *we are not yet accounting the Basis Set Superposition Error (BSSE)*, which will be covered in a later tutorial: :ref:`tutorial-basic-basis-extra`. However, it highlights the impact of even a slight optimization of the basis set. :bdg-info:`background` Manual optimization of the cut-off radii --------------------------------------------------------------- Just relying on :ref:`PAO.EnergyShift` and :ref:`PAO.SplitNorm` is surely convenient when we have to work on a large amount of very different systems. But if our focus is to get the best quality results (for example when working with very subtle effects in energy or fine differences in band structure), then we need to go a step further and optimize the cut-off and matching radii manually. In the following sections, we will start with a basis set generated with a value of 0.001 Ry for :ref:`PAO.EnergyShift` and the default value for the :ref:`PAO.SplitNorm`. However, this is just the starting point, we will be modifying these values the values for the orbital radii manually. By the end of this tutorial, the energy shift will become redundant, and can be actually removed from the input file. But first, we need to have a look at the :ref:`PAO.Basis` block. Now, instead of defining the cardinality just as :ref:`PAO.BasisSize`, we will make use of the :ref:`PAO.Basis` block which expands the information as follows:: %block PAO.Basis H 1 n=1 0 2 P 1 # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals 0.0 0.0 # rcut for first zeta, rmatch for second zeta O 2 # Number of different l n=2 0 2 # n, l (=0, s-orbital), number of zetas 0.0 0.0 # rcut for first zeta, rmatch for second zeta n=2 1 2 P 1 # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals 0.0 0.0 # rcut for first zeta, rmatch for second zeta %endblock PAO.Basis Notice how the H atom has a single 1s level (n=1/l=0) with 2 functions (zetas) and a single set of polarization orbitals (P 1). Meanwhile, Oxygen has two levels, 2s and 2p (n=2/l=0 and n=2/l=1), with an additional set polarization orbitals on the 2p level (P 1). In this example, the input is the same as specifying the flag :ref:`PAO.BasisSize` *DZP*, since we have a double-zeta basis set for all valence orbitals, and we are including a single polarization orbital set for each atom. The values for rcut and rmatch are all set to 0.0; this tells SIESTA that it should automatically find appropriate values for each of those radii using the energy shift. By the end of this tutorial, we will have set each of these radii manually. It is worth noting that all values for radii here should be in **Bohr**. .. note:: As with the first part of the tutorial, in real applications you should converge all rcut values either iteratively or simultaneously. :bdg-success:`basic` Optimizing the first-zeta cut-off radii ------------------------------------------------------------ .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *03-OptimizeFirstZeta* You will now need the following input files: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`water.fdf` .. dropdown:: A peek at the fdf file. .. literalinclude:: work-files/03-OptimizeFirstZeta/water.fdf :linenos: We will first start by optimizing the cut-off radii for the first zeta orbitals. This means we will be changing the first **0.0** we find for each of the orbitals (*the * * *are for emphasis and should not be included in the actual input file*):: %block PAO.Basis H 1 n=1 0 2 P 1 # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals\\ *0.0* 0.0 # rcut for first zeta, rmatch for second zeta O 2 # Number of different l n=2 0 2 # n, l (=0, s-orbital), number of zetas *0.0* 0.0 # rcut for first zeta, rmatch for second zeta n=2 1 2 P 1 # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals *0.0* 0.0 # rcut for first zeta, rmatch for second zeta %endblock PAO.Basis *Modify the value for rcut in the 1s orbital of Hydrogen, going from 4.0 Bohr to 6.0 Bohr by 0.5 Bohr increases. Take note of the basis enthalpy of the system for each value of rcut.* .. dropdown:: Answer: Basis enthalpy for the cutoff radius of Hydrogen 1s :color: success .. code-block:: 4.0 -482.533803 eV 4.5 -482.552463 eV <== 5.0 -482.511975 eV 5.5 -482.480837 eV 6.0 -482.463814 eV Now we choose the cut-off radii for which the basis enthalpy is minimized. Assuming the value is 4.5 Bohr, your basis block will look somewhat similar to this:: %block PAO.Basis H 1 n=1 0 2 P 1 # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals 4.5 0.0 # rcut for first zeta, rmatch for second zeta O 2 # Number of different l n=2 0 2 # n, l (=0, s-orbital), number of zetas 0.0 0.0 # rcut for first zeta, rmatch for second zeta n=2 1 2 P 1 # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals 0.0 0.0 # rcut for first zeta, rmatch for second zeta %endblock PAO.Basis *After choosing a value for the Hydrogen 1s orbital, do the same procedure for the 2s orbital of Oxygen, and then for the 2p orbital. Explore values from 4.5 Bohr to 7.5 Bohr for both orbitals.* .. dropdown:: Answer: Basis enthalpy for the cutoff radius of Oxygen 2s :color: success Having chosen 4.5 for the Hydrogen 1s orbital:: 4.5 -482.554539 eV 5.0 -482.567785 eV 5.5 -482.570432 eV <== 6.0 -482.568240 eV 6.5 -482.563629 eV 7.0 -482.558274 eV 7.5 -482.552222 eV .. dropdown:: Answer: Basis enthalpy for the cutoff radius of Oxygen 2p :color: success Having chosen 5.5 for the 2s orbital:: 4.5 -482.267679 eV 5.0 -482.451362 eV 5.5 -482.551777 eV 6.0 -482.603669 eV 6.5 -482.626368 eV 7.0 -482.631294 eV <== 7.5 -482.625763 eV *Do the values for the optimized cut-off radii coincide with what you would expect from chemical intuition? Which is the largest orbital and which one is the shortest?* .. dropdown:: Answer: The final result :color: success Our PAO.Basis block now looks somewhat like this:: %block PAO.Basis H 1 n=1 0 2 P 1 4.5 0.0 O 2 n=2 0 2 5.5 0.0 n=2 1 2 P 1 7.0 0.0 %endblock PAO.Basis It is consistent with the physics of each atom: the 1s orbital of H is the shortest, while Oxygen has a larger 2s orbital and an even larger 2p orbital. Note that since we are using water as a model system, this basis set might not be the best for other kinds of systems, such as hydrides, in which the electronic density over the atoms is very different. For the purposes of this tutorial, we will stop here the optimization of rcut. In a more general case, it would be better to re-optimize all the orbitals iteratively to guarantee that the rcut values are truly converged. :bdg-success:`basic` Optimizing the second-zeta matching radii -------------------------------------------------------------- .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *04-OptimizeSecondZeta* The input files are essentially the same as in the previous section: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`water.fdf` .. dropdown:: A peek at the fdf file. .. literalinclude:: work-files/04-OptimizeSecondZeta/water.fdf :linenos: We will now expand the previous block to the one we had at the beginning for the DZP basis set, but instead of zeroes in the rcut, we will place our optimized values for the single zeta basis:: %block PAO.Basis H 1 n=1 0 2 P 1 # n, l (=0, s-orbital), number of zetas, include 1 set of polarization orbitals 4.5 *0.0* # rcut for first zeta, rmatch for second zeta O 2 # Number of different l n=2 0 2 # n, l (=0, s-orbital), number of zetas 5.5 *0.0* # rcut for first zeta, rmatch for second zeta n=2 1 2 P 1 # n, l (=1, p-orbital), number of zetas, include 1 set of polarization orbitals 7.0 *0.0* # rcut for first zeta, rmatch for second zeta %endblock PAO.Basis We will do the same as before, but now changing the value for the matching radii (which are still 0.0 in the block). *Using a similar procedure as before, find the value of rmatch for the 1s orbital of Hydrogen that minimizes the basis enthalpy, evaluating it from 1.5 Bohr to 3.0 Bohr. Do the same for Oxygen 2s and then 2p, changing the values from 2.5 Bohr to 4.0 Bohr for each of them.** .. dropdown:: Answer: Basis enthalpy for the matching radius of Hydrogen 1s :color: success .. code-block:: 1.5 -482.666780 eV 2.0 -482.829897 eV <== 2.5 -482.777909 eV 3.0 -482.666422 eV .. dropdown:: Answer: Basis enthalpy for the matching radius of Oxygen 2s :color: success Having chosen 2.0 for the matching radius of the Hydrogen 1s orbital:: 2.5 -482.831605 eV 3.0 -482.832936 eV <== 3.5 -482.816571 eV 4.0 -482.809207 eV .. dropdown:: Answer: Basis enthalpy for the matching radius of Oxygen 2p :color: success Having chosen 3.0 for the matching radius of the Oxygen 2s orbital:: 2.5 -482.815772 eV 3.0 -482.833345 eV <== 3.5 -482.825561 eV 4.0 -482.818233 eV .. dropdown:: Answer: Final PAO.Basis block. :color: success By picking the matching radii that provide the lowest energies, we arrive to the following :ref:`PAO.Basis` block:: %block PAO.Basis H 1 n=1 0 2 P 1 4.5 2.0 O 2 n=2 0 2 5.5 3.0 n=2 1 2 P 1 7.0 3.0 %endblock PAO.Basis Congratulations! You optimized your basis set for a water molecule. In the same way as in the previous section, we will stop here. However, remember that in actual production runs you should verify that all values for the cut-off radii and matching radii are properly converged. :bdg-success:`basic` Checking the result with the water dimer ------------------------------------------------------------- .. dropdown:: **Are you under a local environment?** :color: primary Enter directory *05-OptimizedWaterDimer* It's the water dimer again! We will use the same input files as before, but now we will be using the optimized basis set with the :ref:`PAO.Basis` block. The input files are: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`dimer.fdf` and :download:`monomer.fdf` .. dropdown:: A peek at dimer.fdf .. literalinclude:: work-files/05-OptimizedWaterDimer/dimer.fdf :linenos: Now is the time to put our optimized basis set to test. Remember the binding energy of a water dimer, i.e.: :math:`E_{binding} = E_{dimer} - 2*E_{monomer}` *Calculate the binding energy of the dimer again. How different is our result when compared to the previous ones? Consider both the total energies for each system, and the total binding energy.* .. dropdown:: Answer: Basis enthalpy for different values of split norm :color: success Using our Optimized :ref:`PAO.Basis` block, we get a total energy of **-482.936001 eV** for the monomer and **-966.119097 eV** for the dimer. This gives a binding energy of **-247 meV**. Let's recover our previous results for comparison:: Basis Set | Monomer /eV | Dimer /eV | Binding /meV -------------------|-------------|-------------|------------- Default | -482.136281 | -964.717861 | -445 Rough Optimization | -482.884135 | -966.016972 | -248 Fully Optimized | -482.936001 | -966.119097 | -247 We can see that the binding energy is pretty much the same as the one we got optimizing only roughly; however, note that the total energy of each of the systems is lower when we fully optimize the :ref:`PAO.Basis` block, which may matter if we need finely-tuned results. :bdg-danger:`advanced` Optional: Optimizing the polarization orbitals --------------------------------------------------------------------- .. note:: This step is tipically not necessary, and it may help only in very specific cases. SIESTA already generates good polarization orbitals. .. dropdown:: **Are you under a local environment?** :color: primary Go to the subdirectory *06-PolarizationOrbitals*. One last step in the orbital optimization could be to optimize the polarization orbitals. This is conceptually a bit different from what we have been doing up to now, even if the procedure is somewhat similar. Here, we will not be changing the cut-off radius of the orbital as above; instead, we will be altering the overall shape of the polarization orbital. For the input files: * Pseudos: :download:`O.psml`, :download:`H.psml` * Input files: :download:`water.fdf` .. dropdown:: A peek at water.fdf .. literalinclude:: work-files/06-PolarizationOrbitals/water.fdf :linenos: For that purpose, we will expand the basis block to add one more orbital set per atom:: %block PAO.Basis H *2* n=1 0 2 # n, l (=0, s-orbital), number of zetas *RCUT1 RMATCH1* # rcut for first zeta, rmatch for second zeta n=2 1 1 *Q 0.0* # n, l (=1, p-orbital), number of zetas, charge value for confinement *RCUT1* # rcut for first zeta, rmatch for second zeta O *3* n=2 0 2 # n, l (=0, s-orbital), number of zetas *RCUT2 RMATCH2* # rcut for first zeta, rmatch for second zeta n=2 1 2 # n, l (=1, p-orbital), number of zetas *RCUT3 RMATCH3* # rcut for first zeta, rmatch for second zeta n=3 2 1 *Q 0.0* # n, l (=2, d-orbital), number of zetas, charge value for confinement *RCUT3* # rcut for first zeta %endblock PAO.Basis Note that, in both cases, the polarization orbital belongs to the following principal number (n=3 for Oxygen, n=2 for Hydrogen). In addition, *all polarization orbitals have the same cut-off radius as the orbital they are polarizing*. So, using the results from the previous sections of this tutorial, our starting :ref:`PAO.Basis` block looks like this: .. literalinclude:: work-files/06-PolarizationOrbitals/water.fdf :linenos: :lineno-match: :start-at: %block PAO.Basis :end-at: %endblock One of the recommended methods to generate appropriate polarization orbitals is the charge confinement method (more details in the manual), but suffice to say that this method relies on seeing how the atomic wavefunctions change in the presence of a charge confinement potential, :math:`{ V(r;Q) = Q / r^{2}}`. In order to optimize our polarization orbitals, we just need to change the value of the charge Q. *Modify the charge used in the Hydrogen 2p polarization orbital (the number following Q) from 3.5 to 6.0, and choose a value that minimizes the basis enthalpy. Do the same for the Oxygen 3d polarization orbital, modifying the charge from 5.0 to 8.0.* .. dropdown:: Answer: Basis enthalpy for Hydrogen 2p polarization orbital :color: success .. code-block:: 3.5 -482.523463 eV 4.0 -482.531367 eV <== 5.5 -482.522689 eV 6.0 -482.510356 eV .. dropdown:: Answer: Basis enthalpy for Oxygen 3d polarization orbital :color: success Having chosen 4.0 for Q in the Hydrogen atom: .. code-block:: 5.0 -482.858799 eV 5.5 -482.872189 eV 6.0 -482.877971 eV 6.5 -482.878803 eV <== 7.0 -482.876406 eV 7.5 -482.871858 eV 8.0 -482.865851 eV .. dropdown:: Answer: Final PAO.Basis block. :color: success If we pick the above values for the charge confinement potential charge, then:: %block PAO.Basis H 2 n=1 0 2 4.5 2.0 n=2 1 1 Q 4.0 4.5 O 3 n=2 0 2 5.5 3.0 n=2 1 2 7.0 3.0 n=3 2 1 Q 6.5 7.0 %endblock PAO.Basis *Repeat the calculations for the water dimer using your newly optimized basis set. How do results change?* Use the same input as in the previous section and modify the :ref:`PAO.Basis` block accordingly. .. dropdown:: Answer: Energy of the water dimer with optimized polarization :color: success Let's recover our previous results for comparison:: Basis Set | Monomer /eV | Dimer /eV | Binding /meV -------------------|-------------|-------------|------------- Default | -482.136281 | -964.717861 | -445 Rough Optimization | -482.884135 | -966.016972 | -248 Fully Optimized | -482.936001 | -966.119097 | -247 Optimized Polariz. | -482.980804 | -966.205087 | -243 Note that we have again a slight increase of the binding energy. Again, this differences are small, so we will not have a true final say until we check :ref:`for BSSE`. .. dropdown:: Additional information: decay parameter. :color: info In addition to the charge value, one can also change a decay parameter :math:`{\lambda}`, changing the charge confinement potential to :math:`{ V(r;Q, \lambda) = Q * e^{-\lambda r} / r^{2}}`. This is simply added as an additional number after Q. For example, ``n=2 1 1 Q 3.0 0.4`` has a value of 0.4 for :math:`{\lambda}` and 3.0 for the charge. The decay parameter is not independent of the value of the charge, therefore both values need to be optimized at the same time (or iteratively) for a given orbital. Refer to the manual for more information on the charge confinement potential.