Generating Amorphous Structures

Introduction

Amorphous solid materials become more and more important for a broad spectrum of technologies. In contrast to crystal structures, which are usually well- defined, it is, however, rather difficult to obtain the amorphous structure of a given material in order to perform atomistic simulations. This tutorial gives you some guidelines about generating amorphous structures in Virtual NanoLab (VNL) and ATK using molecular dynamics (MD) simulations at different levels of precision.

Amorphous structures, i.e. structures without any periodic order of the building atoms or molecules, can be found in various kinds of materials, e.g. glass or polymers. Their mechanical and electronic properties often differ quite substantially from the respective crystal. For example, the thermal expansion coefficient of amorphous silica is about one order of magnitude lower than the corresponding value for quartz-crystals.

Consequently, there is a growing interest in the employment in amorphous materials, to tailor various material and device properties in a desired way. In this effort, atomistic simulations can help to understand the behaviour of amorphous strucures at the nano- and microscale.

It is, however, not always straightforward to obtain a suitable amorphous structure of a given material for such simulations. In contrast to the corresponding crystal structure there is no unique unit cell or supercell representation, as the structure inherently bears large degree of randomness. Instead, one should ideally consider an ensemble of realizations of such structures, which is in reality often prohibited by the computational cost of the required calculations.

In this tutorial, you will learn how you can generate amorphous structures that can serve as representative examples and input configurations for further simulations, based on DFT or classical force field methods.

You will learn how to use classical force field molecular dynamics (MD) simulations to obtain a first starting structure of the amorphous material, which can in a second step be refined using semi-empirical or DFT molecular dynamics, when the focus is set on electronic properties.

As there is in general no universally valid recipe to create amorphous structure, the most suitable protocol for a given material depends on a lot of factors, such as the availability of a reliable empirical potential, the desired size of the supercell, the disposable computational resources, etc. This tutorial will therefore only provide a guideline, rather than a rigid protocol, to generate amorphous structures, which might have to be modified according to the particular application at hand.

You will simulate amorphous silicon dioxide, often referred to as fused silica or fused quartz, which is one of the most widely used amorphous materials, e.g. in conventional or specialized glass and also in the semiconductor industry.

You will also learn, how to create interface structures between crystal and amorphous regions.

If you don’t have much experience in MD simulations, you can start with the Molecular dynamics: Basics, which will help you to become familiar with the MD functionality in VNL-ATK.

Note

The simulations in this tutorial use a comparably small number of steps to demonstrate within a reasonable time the principle of generating amorphous structure. For actual production simulations, in order to achieve more relavble amorphous samples with a low amount of structural defects, you should use much longer simulation times. Typically, a melt-quench cycle comprises at least several hundred picoseconds, or even several nanoseconds. Take a look at the References of this tutorial to get an overview over common simulation protocols.

Amorphous Structure Generation with Classical MD Simulations

The most common way to generate amorphous structures is to mimic the experimental protocol of melting a crystal and cooling it down again, by using MD simulations [2], [3]. Although cooling rates accessible in MD simulations can not be chosen as low as the experimental cooling rates, this procedure gives satisfactory results in many cases. The quality of the resulting amorphous structures of course depends much on the ability of the employed force field to describe the material, also further away from equilibrium conditions.

Setting up the geometry

You will start from a crystalline SiO2 system, specifically a cristobalite structure. In contrast to e.g. quartz, the cristobalite unit cell can easily be transformed into an orthogonal cell, which makes the simulation much more convenient. Moreover, the density of cristobalite is very similar to the experimental density of amorphous silica, being 2.2 g/cm3.

Open the Builder and click Add ‣ From Database and search for “sio2”. Select the cristobalite structure and add it to the Stash. To transform it into an orthogonal cell, click on Bulk Tools ‣ Supercell on the right-hand-side panel. Click Conventional and then Transform. For cristobalite this will create a cubic supercell.

../../_images/Builder_cristobalite.png

The resulting cell is yet too small to be a useful container for an amorphous structure, as the periodic boundary conditions would still enforce a certain amount of crystallinity. Therefore, repeat the cell by opening Bulk Tools ‣ Repeat, and apply a 3 x 3 x 3 repetition. This cell, containing 648 atoms, will be the starting structure for the MD simulations.

Setting up the simulation

Send the configuration to the ScriptGenerator. Add a NewCalculator and an Optimize ‣ MolecularDynamics block to the script panel.

Double-click on the NewCalculator block and select the ATK-ForceField calculator with the Pedone_Fe2_2006 potential [1]. Uncheck Print and Save. Close the NewCalculator widget by clicking Ok.

Note

For ATK-versions older than 2017, the ATK-ForceField calculator can be found under the name ATK-Classical.

Now double-click on the MolecularDynamics block. To melt the crystal, you will start with simulations at a very high temperature of 5000 K. As the primary purpose of this initial simulation is to randomize the arrangement of the atoms, no pressure coupling is needed at this stage. Therefore you should select a Type, which gives you an NVT ensemble (take a look at the Molecular dynamics: Basics for a more detailed documentation of MD simulation types in ATK). A good choice is the Langevin thermostat, as it produces a very tight coupling to the heat bath and therefore prevents single atoms from becoming too hot and possibly destabilizing the system.

Set the Steps to 20 000, the Log interval to 5000, and the filename to save the trajectory to SiO2_5000K_low_density_traj.nc. In the Initial Velocity group box, select the Maxwell-Boltzmann type at a temperature of 5000 K.

Enter the same temperature value of 5000 K in the Langevin group box in the Reservoir temperature field. Leave the remaining settings at their default values. Finally, uncheck the Save and Print boxes (as the trajectory is already being saved).

../../_images/molecular_dynamics_sio2_5000K.png

Use the sendto_icon button to send the script to the Editor in order to make some additional changes.

Modifying the script

In some cases, particularly when you attempt to melt comparably small supercells, you may sometimes encounter difficulties to abandon and randomize the crystal structure even at temperatures far beyond the experimental melting point. A possible reason for this behavior, apart from deficiencies of the employed classical potential, is that crystal structures are are additionally stabilized by the periodic boundary conditions, and, generally, artifacts of the periodic boundary increases with decreasing size of the supercell. Another aspect that has to be kept in mind is that the fixed volume of the cell is likely to produce a large system pressure at high temperatures, which can additionally elevate the melting point temperature or in some cases thermodynamically favor another crystal phase over the liquid phase.

To facilitate melting, you can slightly increase the cell volume to artificially lower the density. In the script this is achieved by multiplying the the lengths of the cell vectors by a factor of e.g. 1.1.

Modify the block in which the lattice is defined, as shown in the following script:

# Set up lattice
# Define a scale factor
scale_cell = 1.1
# Scale all 3 cell vectors by this factor
vector_a = [21.498*scale_cell, 0.0, 0.0]*Angstrom
vector_b = [0.0, 21.498*scale_cell, 0.0]*Angstrom
vector_c = [0.0, 0.0, 21.498*scale_cell]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

The exact value of the scaling factor does not play a crucial role, and can be increased as long as you still detect some residual crystallinity after annealing the system at 5000 K.

Send the modified script to the Job Manager and run the simulation.

It will take a couple of minutes to complete the simulation. If you select the produced trajectory file on the LabFloor* and visualize it with the Movie Tool, you will see the atoms moving around in a chaotic fashion, as they frequently cross the periodic boundaries of the cell. To take a better look at the structure, send the last snapshot of the movie to the Builder via the sendto_icon button in the lower right-hand corner of the Movie Tool.

In the Builder you can fold all atoms, that have crossed the periodic boundaries, back into the simulation cell by clicking Bulk Tools ‣ Wrap ‣ Apply.

The final structure may look similar to the example below.

../../_images/Builder_wrap.png

Due to the artificially decreased density the atoms cluster in certain regions, producing an inhomogeneous distribution. However, the purpose of the procedure, i.e. a complete randomization of the previously crystalline structure has been fulfilled very well.

To recover the original density, send this configuration to the ScriptGenerator and repeat exactly the same steps described above (i.e. prepare a script containing a NewCalculator and a MolecularDynamics block with the described settings and send it to the Editor), but now divide the cell vectors by the same scaling factor 1.1 (respectively the factor you have used to lower the density) instead of multiplying them. This procedure will compress the system and recover the original density of the cristobalite structure. Run this simulation at 5000 K for 20 000 steps to equilibrate the silica melt at the new density, send the final configuration from the Movie Tool to the Builder, and Wrap the configuration.

Cooling the structure and adjusting the density

If you inspect the resulting structure with respect to the local atomic arrangements, you will encounter a lot of structural defects, such as undercoordinated atoms (e.g. oxygen atoms with only one silicon neighbor).

This is due to the fact that the resulting amorphous system represents a possible snapshot of the material in the liquid state at very high temperatures. The ultimate goal is usually the generation of solid amorphous material, i.e. at temperatures well below the glass-transition temperature T:sub:`g`. Therefore, you will perform a series of simulations at decreasing temperatures to cool the system to the desired temperature.

Although real amorphous structures always contain a certain amount of structural defects, the goal of such a cooling simulation is to perform the temperature change as slow as possible to allow the system to reduce the amount of defects as much as possible [3]. Unfortunately, the limited accessible simulation time often forces you to use very high cooling rates compared to experimental settings. The settings, suggested in the following, have been chosen to allow for a rather quick demonstration of the principle of the structure generation, which might not be the optimal choice for production simulations. To obtain a more reliable amorphous structure you might consider increasing the simulation times at the individual temperature steps.

Another important parameter is the density of the amorphous structure. There are essentially two different possibilities how to deal with this matter.

One way is to adjust the system volume a priori, to obtain the desired experimental amorphous density and leave the volume at this fixed density throughout the cooling simulations [2]. This technique is the preferred method when the accessible simulation time is very limited, e.g. in DFT-MD simulations. However, it can lead to large residual stress values in the final systems.

The alternative way is to adjust the density continuously during the cooling simulations by coupling the system volume to a barostat at ambient pressure [3]. Here, the resulting density will depend on the quality of the employed potential and you should check, whether it is consistent with experimental values and, in case of large deviations, consider repeating the simulations at fixed volume.

In this example you will use the second technique, i.e. use a barostat to control the system pressure.

Prepare the cooling simulations by sending the last snapshot of the previous simulation to the Builder via the sendto_icon button. Fold the coordinates back into the simulation cell by using the Wrap tool, as described above. Send the configuration to the ScriptGenerator. Add a NewCalculator and an Optimize ‣ MolecularDynamics block. In the NewCalculator widget, use the same settings as in the previous simulations.

Open the MolecularDynamics settings. Select the NPT Melchionna type. Select 5000 steps, a log interval of 500, and set the trajectory name to SiO2_5000K_NPT_traj.nc. Select Maxwell-Boltzmann initial velocities corresponding to 5000 K.

Change the Reservoir temperature to 5000 K, and the Thermostat time scale as well as the Barostat time scale to 20 fs. Set the Bulk modulus to 37 GPa, which corresponds to the bulk modulus of quartz.

Send the script to the Editor.

To cool the system in several steps of decreasing temperature, you can start from the molecular dynamics block, and write a loop over a defined temperature range, as in the follwing example:

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------
# Define a loop over the temperatures
for t in [5000.0, 4000.0, 3000.0, 2000.0, 1500.0, 900.0, 300.0]:
    initial_velocity = MaxwellBoltzmannDistribution(
        # Change the temperature here:
        temperature=t*Kelvin,
        remove_center_of_mass_momentum=True
    )
    method = NPTMelchionna(
        time_step=1*femtoSecond,
        # Change the temperature here:
        reservoir_temperature=t*Kelvin,
        external_stress=1*bar,
        thermostat_timescale=20*femtoSecond,
        barostat_timescale=20*femtoSecond,
        bulk_modulus=370000*bar,
        initial_velocity=initial_velocity,
        mask=[[True, False, False], [False, True, False], [False, False, True]]
    )
    # Define a new trajectory filename for each temperature
    trajectory_filename = 'SiO2_%iK_NPT_traj.nc' % (t)
    md_trajectory = MolecularDynamics(
        bulk_configuration,
        constraints=[],
        trajectory_filename=trajectory_filename,
        steps=5000,
        log_interval=500,
        method=method
    )
    bulk_configuration = md_trajectory.lastImage()

Note that you have to indent all the lines that should be executed within the loop, according to the python language.

Send the modified script to the Job Manager and execute it. It will take some time to complete.

The simulations should run without problems for this system. However, if you encounter extreme, non-stationary cell size fluctuations, or the formation of pronounced vacuum pores during the NPT runs at high temperatures, you can consider switching on the barostat coupling after reaching lower temperatures, e.g. roughly in the region of the suspected glass temperature, while simulating the high temperatures at constant volume.

Analyzing the structure

After the simulations have finished, you can take a brief look at the trajectories using the Movie Tool or the Viewer. You will notice, though, that a lot of atoms have moved beyond the periodic boundaries of the supercell, which makes a visual assessment of the structure difficult. A better visual representation is achieved by sending the final snapshot to the Builder and apply the Wrap tool.

../../_images/Builder_after_cooling.png

If you have a close look, at the local atomic arrangements, you will notice that almost the entire structure is composed of slightly distorted, vertex- sharing tetrahedrons formed by SiO4 units. Additionally, you might encounter some few coordination faults.

For a quantitative assessment, select the SiO2_300K_NPT_traj.nc file on the Lab Floor, and open it with the MD Analyzer tool. Here, you can perform various kinds of analyses of the MD trajectory. For amorphous systems, we are particularly interested in the radial distribution function (RDF). To perform such an analysis, click Load, increase the cutoff radius rmax to 8.0 Å, select Silicon and Oxygen as Element 1, respectively Element 2, and click Plot to display the RDF between these two elements.

Instead of sharp single lines, as in a crystal structure, you will find a distribution of broadened peaks, which is characteristic for an amorphous structure. You can plot the RDF for different element pairs or between all elements.

The amount of structural defects can be investigated by selecting the Coordination number analysis. Select, for example Silicon and Oxygen as Element 1, respectively Element 2, to calculate the coordination of silicon with oxygen atoms. Set the cutoff radius to 2.5 Å, according to the first minima in the RDF curve. Click Plot to create the distribution of coordination numbers. You will see that the majority of silicon atoms possesses its regular coordination number of four, and only a tiny fraction has a lower coordination. Now switch the elements to calculate the coordination number of the oxygen atoms. Click Add Line to add the distribution to the existing plot. Again, the majority of oxygen atoms exhibits its regular coordination number of two, whereas only few atoms are under-coordinated. The MD Analyzer provides a large number of different anaĺysis quantities that can be invoked with a few clicks. Feel free to explore its functions. For example, take a look at the Angular distribution function for Silicon-Oxygen-Silicon angles, or at the Neutron scattering function.

The density of the generated amorphous silica cell can be calculated from the volume of the cell, which can be looked up in the Builder under Bulk Tools ‣ Lattice Parameters, and the numbers of the elements and their respective masses. For the present example, you should obtain a value around 2.6 g/cm3. This value is a little higher than the experimental value of 2.2 g/cm3, which must be considered a small shortcoming of the potential. To obtain exactly the experimental value, you may consider to repeat the cooling simulations under fixed volume, adjusted to reproduce the desired density. The NVT Nose-Hoover Chain thermostat with a linear decrease in temperature might be the right method for such simulations.

The final result can be used as input configuration for subsequent static or MD simulations to calculate further properties of interest. However, in particular with respect to static calculations you should keep in mind, that the generated configuration represents only one possible snapshot out of the amorphous ensemble. To obtain more meaningful results you should consider averaging the calculated values over several different samples.

Refining Amorphous Structures

You can refine the structure obtained via classical MD to compensate some of the possible shortcomings of the force field. This becomes particularly beneficial, if you intend to calculate further properties of the amorphous material based on semi-empirical or quantum mechanical calculations, such as e.g. tight-binding (TB), respectively density-functional-theory (DFT).

The fastest and easiest way is to perform only a TB or DFT geometry optimization to bring the configuration to the nearest energy minimum. This will primarily improve locally distorted geometries, but hardly any possible connectivity defects.

A more profound refinement can be achieved by a short TB/DFT MD annealing and cooling procedure [2]. As these simulation are quite time-consuming they are suited predominantly for smaller systems comprising around 100-200 atoms.

In the following you will be guided through a brief example considering a small silica system with the ATK-DFT calculator. Since the accessible time in DFT-MD simulations is very limited, which makes a reliable constant pressure equilibration very difficult and computationally expensive, you will perform the classical and DFT MD simulations in this particular example at constant volume. The density will remain at the value of the cristobalite crystal, which is very close to the experimental density of amorphous silica.

Setting up the classical simulation

Before running the DFT simulation, you first have to be prepare an amorphous cell using a classical MD simulation. These steps essentially follow the protocol outlined in the previous section.

Prepare a cristobalite cell in the same way as described in the previous section, but use only 2x2x1 repetitions in the Repeat tool of the Builder. This will create a cell containing 96 atoms. Due to this repetition the cell the box size is 14.33 x 14.33 x 7.17 Å3, i.e. the z-direction is thinner compared to the other directions, which might introduce an increased ordering in this direction. To avoid such effects, we will transform the cell to a cubic shape, which has the same volume as the original cell. The final cell size shall therefore be 11.375 x 11.375 x 11.375 Å3. You will perform this shape transformation in two steps during the high temperature simulations where the system is in a liquid state.

The first step comprises setting up a simulation, in exactly the same way as described above, i.e., adding the ATK-ForceField calculator with the Pedone_Fe3_2006 potential and a MolecularDyncamics block with the Langevin type, a step number of 20 000, and a temperature of 5000K. Send the resulting script to the Editor. To achieve a lower density and therefore facilitate the melting, increase the length of the C-vector of the box to its final value of 11.375 Å:

# Set up lattice
vector_a = [14.332, 0.0, 0.0]*Angstrom
vector_b = [0.0, 14.332, 0.0]*Angstrom
# Previous C-vector:
# vector_c = [0.0, 0.0, 7.166]*Angstrom
# New C-vector:
vector_c = [0.0, 0.0, 11.375]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

Run the simulation and send the final snapshot from the Movie Tool to the Builder. Wrap the configuration, set up exactly the same simulation script as in the previous simulation, and send it again to the Editor.

The cell volume has to be reduced to restore the original cristobalite density. Instead of resetting the C-vector to the original value, you can now set the lengths of the A- and B-vectors to 11.375 Å, too, in order to obtain a cubic cell. Run the simulation.

Send the final snapshot to the ScriptGenerator, add the same ATK- Classical claculator as described above, and an Optimize ‣ MolecularDynamics block. In the MolecularDynamics settigns, select the NVT Nose-Hoover Chain type, which allows to use a linear decrease in temperature. Adjust the settings as shown in the following screenshot.

../../_images/molecular_dynamics_nvt_nh.png

These thermostat settings will invoke a linear decrease of the reservoir temperature from 5000 K to 300 K over 50 picoseconds.

Run the simulation. Send the final snapshot via Movie Tool ‣ Send To ‣ Builder ‣ Wrap ‣ Send To to the ScriptGenerator.

Setting up the DFT simulation

In the following you will perform an MD simulation using the ATK-DFT calculator in ATK to anneal the structure and cool it down to room temperature. MD simulations are rather expensive in terms of computational time. In order to generate amorphous structures, it is not absolutely crucial to obtain accurately converged energies. The suggested settings in the present example will therefore employ parameters, whose accuracy is slightly decreased with respect to the default values to accelerate the calculations. However, before you start the DFT simulations, be aware that these simulations can take very long time, up to several days, so it is recommended to run them on a cluster.

To set up the simulation script, add a NewCalculator and two Optimize > MolecularDynamics blocks. In the NewCalculator settings select ATK-DFT calculator. Open the calculator widget to adjust some of the settings. In the Iteration control parameters increase the Tolerance to 0.001 Hartree. In the Basis set/exchange correlation settings, change the exchange correlation Type to GGA. Next, change the basis sets from DoubleZetaPolarized to SingleZetaPolarized.

Depending on your computational resources you may of course choose more accurate settings here, in particular when aiming at production simulations.

Close the NewCalculator settings by clicking Ok.

Add two Optimize ‣ MolecularDynamics blocks to the script. Open the first block. Here, you will set up a short annealing simulation at 3000 K. Therefore, select the Langevin type, set the number of steps to 500 and the log interval to 20. Choose an appropriate trajectory filename, e.g. traj-sio2-DFT-MD-3000K.nc, select Maxwell-Boltzmann initial velocities at 3000 K and a reservoir temperature of 3000 K, as well. Finally, uncheck Save and Print and close the settings.

Open the second MolecularDynamics block. In this simulation you will rapidly quench the system from 3000 K to 300 K. Use again the Langevin Type, set the number of steps to 1000 and the log interval to 20. Select again Maxwell-Boltzmann initial velocities at 3000 K. To cool the system during the simulation set the reservoir temperature to 300 K now, as if the system was transferred from a very hot reservoir to a heat bath at room temperature. Regarding the short simulation time this mehtos is more efficient than using a linear decrease in the reservoir temperature via the NVT Nose-Hoover chain thermostat. Close the widget and run the simulation.

After the simulations have finished, the movie of the second phase, i.e. the quenching simulation, should look similar to the following.

../../_images/MovieTool_DFT_MD.png

The temperature decays from the initial value of 3000 K to room temperature. Of course, this rapid quenching within 1 ps bears a certain deal of approximation compared to the experimental procedure.

Depending on the properties one is interested in calculating, one can continue from the final configuration and perform a DFT geometry optimization.

Creating Crystal/Amorphous Interfaces

A particularly interesting application is given by interfaces between crystal and amorphous structures, as the local properties in proximity of the interface are largely unknown. Such systems can easily be created in the VNL-Builder.

If the crystal is of the same material as the amorphous region, you can start by preparing the crystal configuration. The desired interface plane can be created using the Surface (Cleaver) tool in the Builder. Make a copy of the configuration. The latter system will serve as the basis to create the amorphous structure, as their surface areas already match exactly.

The following example considers the interface between quartz and amorphous silica. The initial structure is built in the Surface (Cleaver) tool by selecting the (10-10) surface and a 4x3 surface cell. Create a Periodic (bulk-like) supercell with a thickness of 12 layers. Perform a geometry optimization of the structure and the cell vectors, using again the ATK-ForceField Pedone_2006_Fe3 calculator, and save the final configuration to a NetCDF file. Drag and drop the file to the Builder, and use it as the initial structure for the amorphous cell.

Prepare the amorphous cell via the same simulation protocol as described above, but apply the pressure coupling only to the Z-axis to preserve the lateral cell dimensions. Unidirectional pressure coupling can be selected in the MolecularDynamics settings, using NPT Melchionna. To couple the barostat only to the Z-component, uncheck all boxes in the Stress coupling group box, apart from the zz box.

After creating the amorphous supercell, you can join the optimized crystal and amorphous cell using the Builders ‣ Interface tool in the Builder panel. Therefore, you can simply drag and drop the configurations from the Stash to First and Second field in the Interface builder. The builder should automatically detect suitable separations of both cells.

../../_images/Builder_interface.png

Click Create to finalize the build process and send the interface configuration to the ScriptGenerator.

Tip

You can learn more about the Interface Builder in the Technical Notes on The Interface Builder in VNL.

If you zoom into the interface, you find encounter a lot of defects, such as dangling bonds. Most of these defects can often be healed by an annealing simulation using classical potentials, even if such potentials are not explicitly designed to include chemical reactions.

Add a NewCalculator, an Optimize ‣ OptimizeGeometry, and an Optimize ‣ MolecularDynamics block. After selecting the desired calculator settings, in this case the Pedone_2006_Fe3-potential, open the OptimizeGeometry widget.

After the interface creation there will probably be a lot of artificially close atomic distances, which will cause a subsequent MD simulation to be unstable. To remove these large force from the interface, you will first run an geometry optimization. The Force tolerance does not need to be extremely low, so you can set a value of 0.5 eV/Å. Increase the Maximum number of steps to 1000. Uncheck the Print box, rename the output file to interface_optimized.nc, and close the OptimizeGeometry settings.

Open the first MolecularDynamics settings. During the MD annealing you should fix the atoms in the central crystal region to prevent the crystal region from melting. To this aim, click on Add constraints to open the Constraints Editor. In the presenter window draw a rectangle around the central atoms of the crystal region to select them, click Add tag from selection, and select Fixed atoms constraints for this group in the table on the side of the widget.

../../_images/constraints_editor_interface.png

Close the Constraints Editor.

In the MolecularDynamics widget, select NVT Nose-Hoover Chain, set the step number to 50 000 and the log interval to 2000 steps. First, you should anneal your structure at a temperature high enough to allow for local rearrangements of atoms, thus facilitating the healing of interfacial coordination defects, but not too high, to prevent melting of the global structure. A temperature of 1500 K is suitable for this purpose. Enter this temperature into the initial velocity temperature and the reservoir temperature fields.

Now add another Optimize ‣ MolecularDynamics block to the script for the cooling simulation. Double-click it and open the Constraints Editor. The previously selected group “Selection 0” should already be present in the table. Select Fixed atom constraints for this group and close the widget. In the MolecularDynamics widget, select the NVT Nose- Hoover Chain thermostat, 50 000 steps, a log interval of 2000 steps, and Configuration velocities as initial velocities. In the thermostat settings set a Reservoir temperature of 1500 K and a Final temperature of 300 K to achieve a temperature decrease dutring the simulation.

Run the simulation by sending the script to the Job Manager.

To finalize the structure, you should run another NPT simulation, this time without any constraints. Use the NPT Melchionna barostat and uncheck all stress coupling boxes except for the zz-box to equilibrate only the cell height perpendicular to the interface plane. Run the simulation for 50 000 steps at 300 K to obtain a room temperature amorphous/crystal interface structure.

Further Examples

Finally, this section provides some further examples of generating amorphpus structres of the technologically relevant materials alumina, titania, and hafnia.

Al2O3

For the creation of an amorphous alumina sample, you can essentially follow the protocol outlined in Refs. [1] and [4]. To use the same classical potential [5], you can either select Matsui_CaMgAlSiO_1994, or, if you use an older version of VNL-ATK where this potential is not available, download the python file attached to this tutorial. If you copy and paste its contents to the top of your simulation script, the potential will be available in the simulation.

To set up the system, you need to import the crystal structure of \(\alpha\)-Al2O3 into the Builder. To this aim, download the Al2O3_corundum_AMS_DATA.cif file. In the Builder, click Add ‣ From Files and open the file. For a more convenient handling in MD simulations the trigonal cell should be converted into an orthogonal cell. The are several ways to achieve this. One possibility is to open the Builders ‣ Surface(Cleave) tool. Select (10-10) Miller indices, click Next and select a 3x1 surface lattice. Click Next and select a Periodic (bulk-like) cell with a thickness of 4. Click Finish to obtain an orthogonal supercell containing 360 atoms.

Now, you should proceed similar to the protocol outlined above, with the exception that the density will be fixed to 3.175 g/cm3 (as proposed in Ref. [4]) during the cooling simulations. The entire protocol reads as follows:

  1. Run a high temperature MD simulation at 5000 K and artificially reduced density.
  2. Cool the system to 3000 K using the NVT Nose-Hoover Chain thermostat.
  3. Compress the system to the desired amorphous density.
  4. Anneal at 3000 K.
  5. Slowly cool the system to room temperature.
  6. Run a final simulation at room temperature to calculate observables.

Send the orthogonal cell to the ScriptGenerator to set up a basic molecular dynamics script. Add a NewCalculator and two Optimize ‣ MolecularDynamics block.

In the NewCalculator block, select ATK-ForceField and the Matsui_CaMgAlSiO_1994 potential.

In the first MolecularDynamics block select the NVT Nose-Hoover Chain thermostat, 150 000 steps, a log interval of 5000 steps and uncheck the Save trajectory box. Set the initial velocities to a Maxwell-Boltzmann- distribution at 5000 K. Set the reservoir temperature to 5000 K and the final temperature to 3000 K. The remaining parameters should be kept at the default values. Finally, uncheck Save and Print, close the calculator widget, and send the script to the Editor.

The second MolecularDynamics block should have similar settings, but 300 000 steps, a reservoir temperature of 3000 K and a final temperature of 300 K.

To obtain an artificially reduced density you should again multiply all three lattice vectors with 1.1, as explained above.

This first MolecularDynamics sequence yields a liquid system at 3000 K and lowered density. Based on this system, the density should be adjusted to achieve the target density of 3.175 g/cm3. As it is often cumbersome to calculate the macroscopic density of an atomic cell and to obtain the exact scaling factors for the cell vectors by hand, this can conveniently be achieved using the ATK unit engine by inserting the following piece of code in the simulation script, right before the second MolecularDynamics block (i.e. after the line bulk_configuration = md_trajectory.lastImage()):

# Get the volume and the lattice vectors of the cell.
V = bulk_configuration.bravaisLattice().unitCellVolume()
lattice_vectors = bulk_configuration.bravaisLattice().primitiveVectors()
# Get the coordinates and the elements.
fractional_coordinates = bulk_configuration.fractionalCoordinates()
elements = bulk_configuration.elements()
# The total and species-resolved number of atoms.
N = len(elements)
N_Al = len(numpy.where(numpy.array(elements)==Aluminum)[0])
N_O  = len(numpy.where(numpy.array(elements)==Oxygen)[0])
# Calculate the current density
density = (N_Al*Aluminium.atomicMass() + N_O*Oxygen.atomicMass())/V
# Get the density in g/cm**3
density = density.inUnitsOf(kiloGram/Meter**3)/1000.0
target_density = 3.175
# Calculate the scaling factor for the lattice vectors.
scale_factor = (density/target_density)**(1.0/3)
# Scale the lattice vectors.
lattice_vectors = lattice_vectors*scale_factor
# Define a new cell.
new_lattice = UnitCell(lattice_vectors[0], lattice_vectors[1], lattice_vectors[2])
# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=new_lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )
# The calculator is still defined and can readily be attached.
bulk_configuration.setCalculator(calculator)

These lines will extract the volume, calculate the density and the scaling factor required to obtain the desired target density of 3.175 g/cm3. Finally, the lattice vectors are scaled accordingly and a new configuration is created, on which the second cooling sequence is performed.

Save the entire script to a python file (an example for the final script can be found in the attachment), and run it via the JobManager or in a terminal using atkpython. It will take a couple of hours to complete. After the simulations have finished you can analyze the structure, as usual.

The complete script can be found in the file script_al2o3_amorphous.py.

TiO2

To simulate amorphous titanium oxide, you will the use the potential of Matsui and Akaogi [6]. This potential has been used in many simulation studies of crystalline and amorphous titania systems, and yields properties in very good agreement with DFT MD simulations. This potential is not yet available as a predefined potential in VNL, but it can easily be set up in a script. Simply download the file Matsui_TiO.py, which is supplied as an attachment to this tutorial. Copy and paste its contents to the top of the script the you will create in the following.

To set up the system open the Builder and use Add ‣ From Database to add the TiO2 rutile crystal structure. Use the Bulk tools ‣ Repeat function to create a 4x4x6 repetition of the system, which contains 576 atoms. Send the configuration to the ScriptGenerator and set up the same basic molecular dynamics script as in the alumina example. Send the final script to the Editor.

Locate the calculator definition as shown below:

potentialSet = Matsui_TiO_1991() # or any other potential class.
calculator = TremoloXCalculator(parameters=potentialSet)

Replace the line starting with potentialSet = ... with the entire content of the attached file Matsui_TiO.py.

When you adjust the volume to obtain the desired density, make sure that you adapt the script using the correct stoichiometry of TiO2, the right element masses, and the right target density. For amorphous titania a density of 3.80 g/cm3, as reported in Ref. [7], can be used.

The complete script can be found in the file script_tio2_amorphous.py.

HfO2

For hafnia you can use the potential of Wang et al. [8], which is available in VNL-ATK as Wang_HfOZr_2012.

To set up the system add the HfO2 crystal structure from the database. Transform the FCC cell into a conventional unit cell, by opening the Bulk tools ‣ Supercell plugin, and clicking Conventional and Transform. Repeat the system 3x3x3 times to obtain a supercell of suitable size.

Set up the same script as in the previous two examples. Make sure to adapt the system and element specific parameters accordingly and use an amorphous density of 7.97 g/cm3 as suggested in [8].

The complete script can be found in the file script_hfo2_amorphous.py.

References

[1](1, 2) A. Pedone, et al: A New Self-Consistent Empirical Interatomic Potential Model for Oxides, Silicates, and Silica-Based Glasses. J. Phys. Chem. B 110, 11780 (2006)
[2](1, 2, 3) E. A. Chagarov and A. C. Kummel: Generation of Realistic Amorphous Al2O3 And ZrO2 Samples By Hybrid Classical and First-Principle Molecular Dynamics Simulations. ECS Transactions 16, 773 (2008)
[3](1, 2, 3) D. J. Cole, M. C. Payne, G. Csányi, S. M. Spearing, and L. Colombi Ciacchi: Development of a classical force field for the oxidized Si surface: Application to hydrophilic wafer bonding. J. Chem. Phys. 127, 204704 (2007)
[4](1, 2) G. Gutierrez and B. Johansson: Molecular dynamics study of structural properties of amorphous Al2O3. Phys. Rev. B 65, 104202 (2002)
[5]M. Matsui: A transferable interatomic potential model for crystals and melts in the system CaO-MgO-Al2O3-SiO2. Miner. Mag. 58A, 571 (1994)
[6]M Matsui and M. Akaogi: Molecular dynamics simulations of the structural and physical properties of the four polymorphs of TiO2. Mol. Sim. 6, 239 (1991)
[7]D. Mergel et al.: Density and refractive index of TiO2 films prepared by reactive evaporation. Thin Solid Films 371, 218 (2000)
[8](1, 2) Y. Wang, F. Zahid, J. Wang, H. Guo: Structure and dielectric properties of amorphous high-k oxides: HfO2 , ZrO2 , and their alloys. Phys. Rev. B 85, 224110 (2012)