2.9. GFN1-xTB: Lithium fluoride

This example refits the Li and F repulsive potentials of the GFN1-xTB model, available in AMS under the DFTB engine. With GFN1-xTB, a lattice optimization of LiF in the rock salt structure collapses the unit cell, because the Li-F interaction is not repulsive enough. This motivates the reparametrization.

../../_images/xTB_LiF_toc_volumescans.png
Table 2.1 In the calculation of the below numbers, the experimental lattice parameters for LiF and Li were used
  Expt. GFN1-xTB Reparametrized xTB (this tutorial)
LiF enthalpy of formation (eV) -6.394 -8.143 -6.393

The fitting is done to better reproduce:

  • DFT-calculated energy-volume curve of LiF
  • Optionally, the stress tensors for each point along the energy-volume curve
  • experimental enthalpy of formation for LiF [reaction enthalpy for ½ F₂(g) + Li(s) → LiF(s)]
  • experimental bond length of F₂ (1.41 Å).

The tutorial will also teach you how you can use different k-space samplings (GFN1-xTB settings) for different jobs.

Note

The input files have already been prepared. If you simply want to run the parametrization without preparing the input, copy the contents of $AMSHOME/scripting/scm/params/examples/xTB_LiF to a new empty directory.

To open the example in the graphical user interface (GUI): SCM → ParAMS, File → Open the job_collection.yaml file in the new directory.

You can then run the parametrization.

2.9.1. Energy-volume scan reference calculation for LiF

Note

The calculation below uses the BAND periodic DFT engine. If you do not have a license for BAND, or want a faster reference calculation for this tutorial, you could instead run the reference calculation using a different DFTB model (DFTBPanel panel, Model → SCC-DFTB, Parameters → QUASINANO2015). You would then fit one DFTB model (xTB) to reference data calculated with another (SCC-DFTB).

First, import the structure and set the DFT settings in AMSinput:

Import LiF_primcell.xyz into AMSinput
Switch to the BANDPanel panel
Task → PES Scan
XC Functional → GGA-D → PBEsolD3(BJ)
Basis set → DZP
../../_images/GFN_amsinput_main.png

Note

For energy-volume scans it is often a good idea to use more dense k-space sampling. Here, we keep the numerical quality at Normal only to have the calculation finish faster. (If you want an even faster calculation, set Numerical Quality to Basic).

Then, set up the energy-volume scan:

Model → Geometry Constraints and PES Scan
Click the AddButton button next to Cell volume
Scale the volume between 0.8 and 1.2 (between 80% and 120% of the original volume)
Set the number of points for SC-1 to 5
../../_images/GFN_amsinput_pesscan.png

Optionally (this will make the calculation more expensive but provide more training data), calculate the stress tensor for each volume:

On the Geometry Constraints and PES Scan panel, check Results for all PES Points.
Properties → Gradients, Stress tensor
Check the stress tensor box.

Save with the name volumescan.ams and run the calculation. Note that it may take about an hour to finish, or more if you enabled the stress tensor calculation (depending on your hardware).

Tip

While the calculation is running you can continue with the step Import experimental formation enthalpy into ParAMS.

When the calculation has finished, see the energy as a function of volume by either

  • opening the calculation in AMSmovie, or
  • selecting the job in AMSjobs, and running Tools → Build Spreadsheet
../../_images/GFN_amsmovie.png

The function run_reference in generate_input_files.py sets up and runs the job via PLAMS:

def run_reference(xyz_file):
    """ 
        Sets up a volume scan of LiF with PBEsol-D3(BJ).
        For each point, the stress tensor is calculated.
        Returns the path to the ams.rkf file

        NOTE: This is an expensive calculation that might take several hours.
    """
    init()
    mol = Molecule(xyz_file)
    s = Settings()
    s.input.band.KSpace.Type = 'Symmetric'
    s.input.band.KSpace.Quality = 'Normal' # N.B. this should be set to 'Good' for production calculations!
    s.input.band.XC.GGA = 'PBEsol'
    s.input.band.XC.Dispersion = 'Grimme3 BJDamp'
    s.input.band.Basis.Type = 'DZP'
    s.input.ams.Task = 'pesscan'
    s.input.ams.PESScan.ScanCoordinate.NPoints = 5
    s.input.ams.PESScan.ScanCoordinate.CellVolumeScalingRange = '0.80 1.20'
    s.input.ams.PESScan.CalcPropertiesAtPESPoints = 'Yes' # needed to get the stress tensor at each point
    s.input.ams.GeometryOptimization.MaxIterations = 0
    s.input.ams.GeometryOptimization.PretendConverged = 'Yes'
    s.input.ams.Properties.StressTensor = 'Yes'
    job = AMSJob(settings=s, molecule=mol, name='volumescan')
    job.run()
    finish()

    return job.results.rkfpath()

2.9.2. Import the energy-volume scan into ParAMS

First, open ParAMS in GFN1-xTB mode:

SCM → ParAMS
Parameters → GFN1-xTB
../../_images/GFN_empty.png

Then, import the results from the reference job:

Job → Add Job
Select the ams.rkf file from the reference job (inside the directory volumescan.results)
This brings up a popup dialog
Use ResultsImporter → Add PES Scan SinglePoints
Task (for new job) → SinglePoint
Properties → Tick relative_energies
If you calculated stress tensor: Properties → Tick stresstensor_diagonal_3d. Note: You can choose several entries in the drop-down list!
../../_images/GFN_add_volumescan.png
Click Import

This creates 5 singlepoint jobs to be run during the optimization.

../../_images/GFN_added_volumescan.png

The training set contains 4 relative energies: volumescan_frame001-volumescan_frame003, …, volumescan_frame005-volumescan_frame003. The energies are taken relative to frame 3, which is the lowest energy frame.

It also (optionally) contains 5 stress tensors, one for each frame. The stresstensor_diagonal_3d property only extracts the diagonal elements of the stress tensor for a 3D system. To also include the off-diagonal elements, you could have checked the stresstensor_3d option in the importer dialog.

Note

It is currently not possible to set a custom unit for the stress tensor in the GUI, but it is possible using scripting. The default unit for the stress tensor is hartree/bohr³.

Save your work before continuing:

File → Save as
Save with the name LiF.params

The function add_volume_scan in generate_input_files.py uses a Results Importer as follows:

def add_volume_scan(results_importer : ResultsImporter, ref_rkf : str):
    """
        The volume scan is added using the add_pesscan_singlepoints() method, which supports extracting the stress tensor.

        If no stress tensor is desired, an alternative is to instead use the add_singlejob() method with task='PESScan'

        See the documentation for more details.
    """
    results_importer.add_pesscan_singlepoints(ref_rkf, properties={
        'relative_energies': {
            'relative_to': 'min',
            'unit': 'eV',
        },
        'stresstensor_diagonal_3d': {
            'unit': 'GPa',
        },
    }, extra_engine='kspace_normal_symmetric')

Here, the extra_engine='kspace_normal_symmetric' option is used to specify job-dependent engine settings.

Note

There are two different params results importers that work on PES scans: add_singlejob and add_pesscan_singlepoints.

The stress tensor can only be imported using add_pesscan_singlepoints and is the results importer used above.

If you did not calculate the stress tensor at each point, you can use the add_singlejob results importer with the pes extractor. That will produce the types of energy-volume curves shown in the ReaxFF: Adsorption on ZnS(110) tutorial. For more information, see PES scans: bond scan, angle scan, or volume scan.

2.9.3. Import experimental formation enthalpy into ParAMS

ParAMS can be used to train to experimental reference data. In that way, you do not need to run expensive DFT reference calculations. Here, we demonstrate it for LiF. The reaction

½ F₂ (g) + Li(s) → LiF(s)

has ΔH⁰ = -617 kJ/mol = -147 kcal/mol = -0.2350 Ha = -6.394 eV. Let’s add this reaction energy to the training set.

First, let’s add a geometry optimization job for F₂(g), and two single point jobs for Li(s) and LiF(s). Li has the body-centered cubic (bcc) structure, and LiF has the NaCl (rock salt) structure.

Open AMSinput, and switch to the DFTBPanel panel.
Task → Geometry Optimization
Click the MoreBtn next to Geometry Optimization
Max iterations → 30
../../_images/GFN_geoopt_settings.png
Details → Expert AMS
Scroll down to GeometryOptimization, and tick the Pretend converged box
../../_images/GFN_geoopt_expertams.png
Draw an F₂ molecule: in the 3D molecule area: Type f and click the positions of the two atoms. Or use the XTool periodic table to select fluorine.
Click the Pre-optimize button on the bottom right of the main panel, and wait for the optimization to finish
../../_images/GFN_F2_preopt.png
File → Add to ParAMS (Ctrl-T). This brings up the ParAMS window.
In the ParAMS window, go to the Jobs panel and change the JobID of the added job to F2.
../../_images/GFN_rename_F2.png

The added job F2 is a Geometry Optimization job, meaning that the geometry will be optimized for every set of parameter values during the parametrization.

We now perform similar steps for Li and LiF, but add them as SinglePoint jobs:

Go back to AMSinput
Delete the F₂ molecule
On the main panel, set Task → SinglePoint
Edit → Crystal → Cubic → bcc
In the Presets dropdown, select Li and click OK
../../_images/GFN_preset_Li.png
File → Add to ParAMS (Ctrl-T). This brings up the ParAMS window.
In the ParAMS window, change the JobID of the added job to Li_bulk.
Go back to AMSinput
Delete the Li crystal
Edit → Crystal → Cubic → NaCl
In the Presets dropdown, select LiF and click OK
../../_images/GFN_preset_LiF.png
If you get a question about mapping to unit cell, select No
File → Add to ParAMS (Ctrl-T). This brings up the ParAMS window.
In the ParAMS window, change the JobID of the added job to LiF_bulk.
Go back to AMSinput and close the AMSinput window
Go back to the ParAMS window
../../_images/GFN_added_LiF.png

Note that the Reference Engine column says DFTB - it corresponds to the settings used in AMSinput. It can be used to use ParAMS to calculate reference values. However, in this tutorial we will use experimental reference data. Therefore, let’s delete the reference engine from the three new entries.

Double-click in the Details for the F2 job (where it says Geometry Optimization + ...).
In the Reference Engine ID drop-down, choose None
../../_images/GFN_unset_ref_engine.png
Set the Reference Engine ID to None also for the Li_bulk and LiF_bulk jobs.
../../_images/GFN_unset_all_ref_engines.png

Now you can add a reaction energy:

In the ParAMS window, switch to the Training Set panel
In the menu, choose Training Set → Add → Energy (Ctrl-E).
../../_images/GFN_added_empty_energy.png
Double-click the new Energy row in the Detail column.
In the Energy box, type 1.0*energy("LiF_bulk")-energy("Li_bulk")-energy("F2")
../../_images/GFN_prebalance.png

Tip

When typing energy expressions, you can use tab-completion. For example, if you start typing energy("Li and press TAB multiple times, you will cycle through energy("Li_bulk") and energy("LiF_bulk").

Click the Balance button. This balances the chemical reaction to match the given coefficient of 1.0.
In the Values box, type the reaction energy. The needed unit is shown above the box (the unit is controlled by your preference). The value is -617 kJ/mol = -147 kcal/mol = -0.2350 Ha = -6.394 eV
../../_images/GFN_postbalance.png
Click OK
../../_images/GFN_added_enthalpy.png
Save your work before continuing.

The function add_enthalpy_of_formation in generate_input_files.py uses the Results Importer add_reaction_energy. The needed .xyz files can be found in $AMSHOME/scripting/scm/params/examples/xTB_LiF.

def add_enthalpy_of_formation(results_importer : ResultsImporter) :
    """
        Adds experimental enthalpy of formation of LiF
        The jobs are added with Good k-space quality, since total energies of different
        supercells are sensitive to the k-space sampling. (Moreover Li is a metal)
    """
    LiF = Molecule('LiF_primcell.xyz')
    Li = Molecule('Li_primcell.xyz')
    F2 = Molecule('F2.xyz')
    results_importer.add_reaction_energy(
        reactants=[Li, F2],
        reactants_names=['Li_bulk', 'F2'],
        products=[LiF],
        products_names=['LiF_bulk'],
        normalization='p0', # reaction energy per LiF (first product)
        task='GeometryOptimization',
        reference=-6.394,
        unit='eV',
        extra_engine='kspace_good_symmetric' # applied to the Li_bulk and LiF_bulk jobs. The F2 molecule is nonperiodic so the k-space settings are ignored
    )

The extra_engine option specifies the DFTB settings to be used during the parametrization.

Important

When adding the reference value manually, verify that the stoichiometric coefficients are correct. Here, the Li_bulk job contains one formula unit of Li (the primitive Li unit cell), the LiF_bulk job has one formula unit of LiF (the primitive LiF unit cell), and F2 has one formula unit of F₂. This means that the chemical systems perfectly match the species in the chemical equation.

However, if you for example were to import the conventional unit cell of LiF, which has 4 LiF formula units (Li₄F₄), then the reaction would be written as

½ F₂(g) + Li(s) → ¼ Li₄F₄(s), ΔH⁰ = -617 kJ/mol

or

2 F₂(g) + 4 Li(s) → Li₄F₄(s), ΔH⁰ = -2468 kJ/mol

Note

The above box shows two different ways of writing the same reaction, with different reaction energies. During the parametrization, the error of the predicted reaction energy will be 4 times larger if you use the second equation (ΔH⁰ = -2468 kJ/mol). That will implicitly make this training set entry more important, at the expense of other training set entries.

Set Sigma to what you think would be an acceptable error for your reaction energies, and set Weight to how important you consider the training set entry to be.

See also: Sigma vs. weight: What is the difference?.

2.9.4. Import experimental F₂ bond length

For the F2 geometry optimization job from the previous section, we can add another training set entry.

In the ParAMS window, on the All panel, select the F2 job.
Training Set → Add → Bonds
On the Training Set panel, double-click in the Details column of the newly added entry
Change the value to 1.41 (angstrom).
Click OK.
Save your work before continuing.
../../_images/GFN_added_bond.png

The add_bond_length function in generate_input_files.py:

def add_bond_length(results_importer : ResultsImporter):
    """ add experimental F2 bond length. This requires that the job F2 has been defined (which is done in add_enthalpy_formation) """
    results_importer.data_set.add_entry('distance("F2",0,1)', reference=1.41)

This function adds an entry to the dataset with the distance extractor. See also: Demonstration: Working with a DataSet

2.9.5. Job-dependent engine settings (k-space)

GFN1-xTB is an electronic structure method. For periodic crystals, it is important to converge the k-space sampling, especially if the energies of the jobs are used. In this tutorial, both the energy-volume scan and the formation enthalpy make use of the energies, and compare energies for different lattices. What should the k-space samplings be for those jobs?

It is up to you to decide. Here, we will demonstrate how to set a Good k-space quality for the Li_bulk and LiF_bulk jobs, and how to use Normal k-space quality for the volume scan jobs. For F2 (gas-phase), the k-space settings are ignored.

We set KSpace Type to Symmetric, because the systems (primitive unit cells of Li and LiF) are very symmetric (it is faster than a Regular grid).

See also

k-space documentation for more details about k-point grids (Regular, Symmetric, Quality, NumberOfPoints, …)

In the GUI, the engine settings during the parametrization are referred to as “ParAMS Engines” (in contrast to settings used for reference calculations, that are referred to as reference engines).

There is always a ParAMS engine called ParAMS, which uses the default engine settings.

We will create two new ParAMS engines

  • kspace_normal_symmetric : KSpace Type Symmetric and KSpace Quality Normal
  • kspace_good_symmetric : KSpace Type Symmetric and KSpace Quality Good
In the ParAMS window, go the the Jobs panel
The ParAMS Engine column has the value ParAMS for all jobs.
Select the 5 volumescan_frame00* jobs. First select the top one, then hold Shift and select the last one.
../../_images/GFN_select_volumescans.png
Double-click in the Detail column
In the ParAMS Engine ID dropdown, select New ParAMS Engine. This will create a new engine called ID-1.
../../_images/GFN_new_params_engine.png
Click the AMSinput button to the right
This brings up a warning stating that only the engine settings will be used. Click OK.
On the DFTBPanel main panel, set Periodicity to Bulk. This activates the k-space options.
Click the MoreBtn arrow next to K-space.
Set K-space grid type to Symmetric
../../_images/GFN_kspace_symmetric.png
Close the AMSinput window (File → Close Job)
Pass to ParAMS? Yes
Click OK in the Job Details window

In the Jobs panel, the ParAMS engine is now given as ID-1 for the volumescan jobs.

../../_images/GFN_added_id1.png

We can give the engine a more descriptive name:

Go to the Engines panel
../../_images/GFN_before_rename_id1.png
Double-click on ID-1
Change the Engine ID from ID-1 to kspace_normal_symmetric
../../_images/GFN_rename_id1.png
Click OK

On the Jobs panel, the ParAMS Engine ID has also been updated:

../../_images/GFN_renamed_id1_jobs.png

Warning

When you select many jobs, double-click in the details column, and make any changes, all settings (AMS settings, Reference Engine ID, and ParAMS Engine ID) will be copied to the selected jobs.

Tip

If you have many jobs, use the search box at the bottom of the table to filter. For example, type volumescan to only see the jobs matching volumescan. Remember to clear the search box to see all jobs again.

Repeat the above steps for the LiF_bulk and Li_bulk jobs. Select the jobs, double-click in the Details column, create a new ParAMS engine, set K-space grid type to Symmetric and the Quality to Good. Name the engine kspace_good_symmetric.
../../_images/GFN_kspace_good.png ../../_images/GFN_final_job_setup.png
Save your work before continuing.

The function add_engines in generate_input_files.py adds two new engines to the Engine Collection:

def add_engines(results_importer : ResultsImporter):
    """ Control the k-space settings for each job by adding different engine definitions. """
    s = Settings()
    s.input.dftb.KSpace.Type = 'Symmetric'
    s.input.dftb.KSpace.Quality = 'Normal'
    results_importer.job_collection.engines.add_entry('kspace_normal_symmetric', Engine(s))

    s = Settings()
    s.input.dftb.KSpace.Type = 'Symmetric'
    s.input.dftb.KSpace.Quality = 'Good'
    results_importer.job_collection.engines.add_entry('kspace_good_symmetric', Engine(s))

The labels 'kspace_normal_symmetric' and 'kspace_good_symmetric' can be referred to using the extra_engine argument when calling a results importer. See the previous examples above.

2.9.6. Set parameters to optimize and their ranges

In xTB, the repulsive potential is controlled by the two parameters alpha and Z for each element.

At the bottom of the Parameters panel, type F: in the search box.
Tick the Active checkbox for F:Z and F:alpha.
../../_images/GFN_activate_F_parameters.png
Type Li: in the search box
Tick the Active checkbox for Li:Z and Li:alpha.
Clear the search box

Optionally, you can also change the minimum and maximum allowed values of the parameters by changing the numbers in the Min and Max columns.

Double-check that only the four parameters F:Z, F:alpha, Li:Z, and Li:alpha are active:

In the Active dropdown list, untick Show Non-Active.
This should show only the 4 parameters (make sure the search box is empty).
../../_images/GFN_active_parameters.png

The function generate_parameter_interface in generate_input_files.py creates the parameter interface.

def generate_parameter_interface():
    """ 
        Generate parameter_interface.yaml activating 
        the Li and F repulsive potentials for optimization 
    """
    # s is the *default* engine settings used during the parametrization
    # It is normal to leave it as an empty Settings and set the engine settings per job
    s = Settings() 
    interface = GFN1xTBParameters(settings=s)

    # Repulsive potential of Li.
    interface['Li:alpha'].is_active = True
    interface['Li:Z'].is_active = True

    # Repulsive potential of F.
    interface['F:alpha'].is_active = True
    interface['F:Z'].is_active = True

    interface.yaml_store("parameter_interface.yaml")

2.9.7. Set the optimizer settings

Switch to the Settings panel
In the Settings dropdown, choose Optimizer: CMAOptimizer
Set popsize to 8
Set sigma to 0.05
In the Settings dropdown, choose Logger
Set logger_every to 5
In the Settings dropdown, choose Max Evaluations
Set max_evaluations to 400
../../_images/GFN_settings.png
Save before continuing

The params.conf.py file contains:

optimizer = CMAOptimizer(popsize=8, sigma=0.05)
logger_every = 5
max_evaluations = 400

2.9.8. Run the xTB parametrization

In the ParAMS window, select File → Run to run on the local machine

Alternatively, select the job in AMSjobs and first choose which queue to run on, before submitting the job.

Use the following command in the directory containing job_collection.yaml, job_collection_engines.yaml, parameter_interface.yaml, training_set.yaml, and params.conf.py:

"$AMSBIN/params" run

Depending on your hardware, after about 30-60 minutes you should have obtained a good set of parameters.

Tip

The results will automatically update in the GUI while the parametrization is running.

If you think the results are “good enough”, you do not need to wait for the parametrization to finish. In the GUI, select the job in AMSjobs and choose Job → Kill to stop the job.

2.9.9. Results of the xTB parametrization

For example,

  • The optimized parameters are shown on the Parameters tab.
  • On the Graphs tab you can see the running loss function (with the CMA-ES optimizer running in parallel, some of the evaluations are logged out of order - that is completely normal):
../../_images/GFN_running_loss.png
  • On the Graphs tab you can create a scatter plot of predicted vs. reference stress tensor components, here compared to the original GFN1-xTB:
../../_images/GFN_prediction_stress.png

Tip

You double-click on a plot axis to access plot configuration settings.

  • On the Results tab you can select Training best: energy to see the energies
../../_images/GFN_predictions_energy.png
  • The xtb_files directory contains all the files needed to run new AMS calculations with the optimized parameters. On the DFTBPanel panel in AMSinput, set Model to GFN1-xTB and Parameters directory to the xtb_files directory. If scripting, set Engine DFTB%ResourcesDir to the xtb_files directory.
  • xtb_files/elements.xtbpar contains the optimized values of the alpha and Z parameters:
#Element eta Gamma alpha Z k_XB Electronegativity
H 0.470099 0.0 2.2097 1.116244 0.0 2.2
He 1.441379 1.5 1.382907 0.440231 0.0 3.0
Li 0.205342 1.02737 0.3471898244638939 2.481445788775993 0.0 0.98
Be 0.274022 0.900554 0.865377 4.07683 0.0 1.57
B 0.34053 1.3 1.093544 4.458376 0.0 2.04
C 0.479988 1.053856 1.281954 4.428763 0.0 2.55
N 0.476106 0.042507 1.727773 5.498808 0.0 3.04
O 0.583349 -0.005102 2.004253 5.171786 0.0 3.44
F 0.788194 1.615037 2.465746060529967 4.812536210787959 0.0 3.98

See also

Getting Started: Lennard-Jones Potential for Argon for more detailed descriptions about the output files and plots.

Note

When using the Add PESScan SinglePoints option, the GUI does not plot any energy-volume curves (and the pes_predictions folder in the output is empty). To get such plots, instead use Add SingleJob with the pes property (see ReaxFF: Gaseous H₂O)

You can find some example output files in $AMSHOME/scripting/scm/examples/xTB_LiF/example_output/training_set_results/best.