In many calculations you can save a lot of time by using a converged calculation as a starting guess for another calculation. In fact, for finite bias calculations of complicated transport system it can sometimes be the only way to obtain stable convergence.

To do this requires a little bit of coding, but normally it is only required to add a few lines to your existing script.

Background

The default start guess for a self-consistent calculation in ATK is based on "neutral atoms", i.e. the basis orbitals are populated as if each atom was isolated in space. This can in many cases be quite far from the "real" state of the atoms, even in a simple crystal, so it seems obvious that you could save a lot of time by providing a better start guess. A typical scenario in which you very simply have access to a better start guess is if you just before performed a calculation for a similar system - this is the case of looping over a particular variable to see how much the results change. Therefore, you can apply this method for instance

• to obtain good convergence at finite bias; the recommended approach to obtain an I-V curve is to "climb" in bias, starting at zero bias and then initialize each finite-bias calculation from the previous, lower bias;
• when checking convergence in various parameters, like k-points or mesh cut-off;
• while scanning a lattice constant for a bulk system, or a coordinate for one or several atoms;
• as a special case: to initialize the anti-parallel configuration in a magnetic tunnel junction.

Since the initial guess is applied to the density matrix, it is a requirement that the original configuration/calculation (to be used to initialize the new one) and the new one have exactly the same sets of basis orbitals. Thus, the two configuration must have

• the same atoms, element by element (the order must not change either), and
• the same basis set, atom by atom.

The atoms do not have to be in the same geometric positions, however. But, if the atoms are moved too far, the initial guess is likely worse than using neutral atoms.

The way to initialize one calculation from another is to add a keyword initial_state to the setCalculator() method on a configuration. The initial state object is typically read from a NetCDF file, or a variable in the script if you are looping over e.g. k-points or bias.

It is probably easiest to understand how to do this by looking at a couple of examples.

Examples

Simple case

For a very simple example, assume you have performed a calculation for a water molecule with a mesh cut-off which is lower than default. Now you want to use this as the initial guess for a higher mesh cut-off, to speed up that calculation. It will be assumed that the first calculation produced a NetCDF file called water_low_accuracy.nc, and that it contains only one MoleculeConfiguration. Then, assume you are now in the process of setting up new calculation, with the higher accuracy, in the Script Generator in VNL. You already have the molecular configuration and the New Calculator blocks defined, as in the figure below.

Now, double-click Initial State in the "Blocks" panel to insert such a block, and then double-click the inserted block in the "Script" panel to open it. Tick the checkbox "Use old calculation", and then use the "..." icon under "Old calculation" to locate the file water_low_accuracy.nc.

The effect of this is precisely to insert an initial_state variable for setCalculator(), as mentioned above:

If you actually do this in practice, for water, using 40 H as the "low accuracy" mesh cut-off and 120 H for "high accuracy", you will find that the first calculation takes 20 seconds (on a particular machine) and runs for 9 self-consistent iterations. The second one (which uses the first one as initial state) converges in only 4 steps and takes 60 seconds. If you don't use the initial state, the high-accuracy calculation also needs 9 iterations, which take 100 seconds.

So,you can in this simple example save about 20% of the calculation time, although it should be noted that the first calculation also takes some time, and there is only an actually reduction of the total calculation time if the starting guess it produces is so good that we reduce the number of iterations with the more expensive settings significantly. Your mileage with this approach (to first do a quick-and-dirty calculation to use as starting guess) will thus vary strongly from case to case.

Stepping up in bias

For the case of an I-V curve (climbing in bias), see the I-V mini-tutorial.

Looping over parameters

In calculations of the type ATK performs, there are numerical parameters that you don't always know a priori how to choose optimally (meaning, where a good balance is found between calculation time and accuracy). This applies to parameters like k-point sampling and mesh cut-off. The only way to find out, besides using some reasonable physics logic to at least get a rough idea, is to try different values and see how the results change. This can often be done on a small test system, and the results can then be applied to the larger, real system of interest.

The script below checks how many k-points are needed to get a good value for the total energy of a gold crystal. To get reliable results when doing parameters scan, it's a good idea to lower the tolerance of the self-consistent loop (increase the accuracy, that is), so we set it two orders of magnitude lower than default in this script.

# Try k-point sampling grids in the range 1x1x1 to 19x19x19
k_values = range(1,20) 

# Create a list to store the computed total energies
energies = []

# The configuration to compute is bulk gold
lattice = FaceCenteredCubic(4.07825*Angstrom) 
elements = [Gold] 
cartesian_coordinates = [[ 0.,  0.,  0.]]*Angstrom
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice, 
    elements=elements, 
    cartesian_coordinates=cartesian_coordinates
)

# For the first step there is no initial state
initial_state = None

# Loop over the different k-point sampling grids
for Nk in k_values:     
    numerical_accuracy_parameters = NumericalAccuracyParameters(
        k_point_sampling=(Nk, Nk, Nk)
    ) 
    iteration_control_parameters = IterationControlParameters(
        tolerance=4e-07,
    )
    calculator = LCAOCalculator(
    numerical_accuracy_parameters=numerical_accuracy_parameters,
        iteration_control_parameters=iteration_control_parameters,
    )
    # This is where we use the previous calculation as starting guess!
    # Note that None is a valid value for the initial_state
    bulk_configuration.setCalculator(         
        calculator,         
        initial_state=initial_state
    )

    # Compute and store the total energy
    energies.append(TotalEnergy(bulk_configuration).evaluate())

    # Save the converged result to use as initial state for the next value of Nk
    initial_state = bulk_configuration

# Print results (on the master node only)
if processIsMaster():
    for Nk,energy in zip(k_values,energies):
        print Nk, energy

Running this script, you will see that a sampling of 12 or 13 k-points along each axis is sufficient for good total energies for gold (for this basis set, etc). If you compare to the same calculation without any initial state, i.e. if each calculation is started from neutral atoms, the number of self-consistent iterations is typically 2 times larger, at least when the number of k-points is large (i.e. the more time-consuming parts of the script), thus the benefit of using an initial state in this type of loop can be quite substantial, and is almost always a good idea.

Another very common task is to loop over a variable, like a lattice constant, to find the lowest total energy for a configuration. The next script will find the equillibrium lattice constant of gold:

# Try lattice constants in steps of 0.05 Angstrom
a_values = numpy.arange(3.9,4.4,0.05)

# Create a list to store the computed total energies
energies = []

# For the first step there is no initial state
initial_state = None

for a in a_values:
    # For each lattice constant, create the configuration
    lattice = FaceCenteredCubic(a*Angstrom)
    elements = [Gold]
    cartesian_coordinates = [[ 0.,  0.,  0.]]*Angstrom
    bulk_configuration = BulkConfiguration(
        bravais_lattice=lattice,
        elements=elements,
        cartesian_coordinates=cartesian_coordinates
    )
    # 13x13x13 is a sufficient k-point sampling for gold
    numerical_accuracy_parameters = NumericalAccuracyParameters(
        k_point_sampling=(13, 13, 13)
    ) 
    iteration_control_parameters = IterationControlParameters(
        tolerance=4e-07
    )
    calculator = LCAOCalculator(
        numerical_accuracy_parameters=numerical_accuracy_parameters,
        iteration_control_parameters=iteration_control_parameters,
    )
    # Initial state is set here
    # Note that None is a valid value for the initial_state
    bulk_configuration.setCalculator(
        calculator,
        initial_state=initial_state,
    )

    # Compute and store the total energy
    energies.append(TotalEnergy(bulk_configuration).evaluate())

    # Save the converged result in a separate variable, to use as initial state for the next value of a
    initial_state = bulk_configuration

# Print results (on the master node only)
if processIsMaster():
    for a,energy in zip(a_values,energies):
        print a, energy

    # Simplistic approach to finding the minimum: fit a parabola
    coefs = numpy.polyfit(numpy.array(a_values),numpy.array(energies),2)
    print "Fitted minimum: a=%s Angstrom" % (-coefs[1]/coefs[0]/2)

Notes

• To optimize the unit cell size for a system like this, it is often more appropriate to use the OptimizeGeometry method in ATK, and let it optimize the cell as well. The calculation above is still useful, however, since it also gives an idea of the error we make by having a non-optimized cell, and the values can furthermore be used to compute the bulk modulus if you do a more proper 3rd degree polynomial fitting!
• In these examples the results are not saved, only printed to the screen/log file. If you want to do some further processing, like compute the band structure for each configuration and compare those, add such statements inside the loop and store the results in a NetCDF file.
• The obtained result is not universal. First of all it's specific for gold, but more importantly you typically need a higher k-point sampling to get accurate forces or stress, even for gold.
• The result of fitting the parabola is a=4.17 Å; if we use GGA instead the value changes to 4.19 Å, i.e. in both cases about 2-3% above the experimental value 4.08 Å. It is not expected that pseudopotential DFT always gives the exact experimental value, so this is not an error. A more relevant comparison is for instance what a corresponding plane-wave calculation gives, and these typically also overestimate the lattice constant of gold by 1-3% (see e.g. J. Am. Chem. Soc. 128, 3659 (2006)).