HRSMC short rotation: Simulations of explosive sensing materials

HRSMC short rotation report

--- ---

Simulations of explosive sensing materials

--- ---

Robert Kamphof

11328649

Dr. D. Dubbeldam

Dr. S. Ingemann Jorgensen

2

Table of Contents

Introduction ... 3

Project overview and methods ... 3

Results and discussion ... 4

Can we trust simulations to be true to reality? ... 5

Errors encountered during the simulation... 5

Structural correlations to predicted pXRD patterns ... 6

Conclusion ... 8

Evaluation and reflection ... 8

Appendix A: optimised structures ... 9

3

Introduction

Metal-organic frameworks (MOFs) are an emerging class of materials that find applications in catalysis, separation and sensing. Consisting of metal ion nodes and organic ‘linker’ molecules, MOFs form porous three-dimensional structures. The choice of metal and linker as well as the synthesis conditions all determine which of a nigh-infinite variety of structures is formed. During my bachelor project, I synthesised a new structure, RK26, for the purpose of sensing TATP, a dangerous homemade explosive that has been used during terrorist attacks such as the attack on Zaventem Airport in 2016. RK26, with molecular formula [Tb(bpdc)1.5(H2O)(DMF)]•2H2O•2DMF, seems to be a

promising candidate for continued research, but being a new structure, much is still unknown about its properties and reactivity. For example, it exhibited significant changes in its pXRD pattern upon contact with different solvents. Whether this indicates poor stability or great structural flexibility is yet unknown and hard to determine. In order to shed light on the structural properties of this material, computational studies could be of great value. Using computational simulations, structures and properties can be predicted

and then compared with experimental data to validate hypotheses. This would be hard to do or even impossible without computational methods. For example, the predicted pXRD pattern of a simulated structure with empty pores can be compared to experimental data in order to determine if the structure has been successfully emptied of its solvent. During this project, I performed quantum optimisations on the structure of RK26 as it was measured via X-ray crystallography in order to find out more about this material and to learn more about computational research in general.

Project overview and methods

All computations were performed using the VASP software package, the Vienna Ab Initio Simulation Package for atomic scale materials modelling from first principles. It was designed to be a simulation environment which is adaptable and facile to use. It specialises in simulations on solid materials such as zeolites, carbon nanotubes and MOFs, which are generally difficult to simulate.

Using VASP, several simulations were performed with the following purposes in mind:

 To obtain the minimum-energy structures of RK26 with filled pores, with non-coordinated solvents removed, and finally the evacuated structure.

 To see the effect of the pore loading on the angles of the unit cell.

 To obtain predicted pXRD patterns of the resulting structures, which can be compared to experimental data, and to gain qualitative insight into the origin of the peaks of the patterns. The starting structure for the simulations was built based on X-ray crystallographic data obtained by the Berkeley X-ray crystallography facility for single-crystal XRD. This structure was modified in the following ways using a structure editing program called Chimera:

 Some molecules were edited to remove mistakes in their molecular structure.  Hydrogen atoms were added where needed.

 Three substructures were then produced:

o One was used as-is with full solvent loading (1)

o One was emptied of all non-coordinated solvents, but coordinated molecules were retained (2)

o Finally, one substructure was completely stripped of all solvent molecules. (3)

Figure 1: The structure of RK26 resulting from crystallography experiments. For clarity, non-coordinated solvents have been removed.

4 The following jobs were then performed:

# Job name Description Starting

structure

Result

1 RK26ISIF0 Experimental structure; cell angles and volume locked

Experimental Successful 2 RK26solventfreeISIF0 Like (1), with uncoordinated

solvents removed

Experimental Successful 3 RK26cleanISIF0 Like (1), with all solvents removed Experimental Successful 4 RK26ISIF3 Like (1), with cell angles and

volume released

(1) result Failed 5 Rk26solventfreeISIF3 Like (2), with cell angles and

volume released

(2) result Failed 6 RK26cleanISIF3 Like (3), with cell angles and

volume released

(3) result Failed 7 RK26ISIF4 Like (1), with only cell angles

released

(1) result Successful 8 RK26solventfreeISIF4 Like (2), with only cell angles

released

(2) result Successful 9 RK26cleanISIF4 Like (3), with only cell angles

released

(3) result Successful 10 RK26anglesto90 Like (1), with cell angles locked to

90 degrees

(1) result Successful 11 RK26cleananglesto90 Like (3), with cell angles locked to

90 degrees

(3) result Successful 12 RK26cleananglesreleased Follow-up of (11) where cell angles

are released

(11) result Successful

Results and discussion

Images of the structures resulting from the simulations can be found in the appendix. A general trend that was observed was that without solvents, linker molecules tended to stick together in tight pairs, while they were stacked more spaciously when solvent molecules were present in the pores. This is easily explained by the solvent molecules interacting with the linkers, weakening the π-π stacking interaction between the linker pairs.

Another trend was that in structures where the coordination sphere of the metals was disturbed (e.g. the structures where

coordinated solvents were removed), the distance between two metal atoms decreased. This could be due to the bond order between the metals increasing as a result of the loss of the solvent molecules, or due to the fact that with those molecules out of the way, there was more room for the metals to move closer. An exception to this trend could be

seen in job (11), where the cell angels were set to 90. During that simulation, the metal atoms for the ‘clean’ structure were pulled apart to the point where Tb-O bonds between the metal and the linker started to break. However, it probable that this observation is due to the unphysical nature of the simulation; because the cell volume and angles were not allowed to change, there was strain on the structure that could not be resolved. As soon as the cell angles were released in job (12), the

Figure 2: In simulations without solvent molecules, linkers tended to stick more closely together.

5 metal atoms snapped together again within a few cycles. One thing to note is that the result of job (12), where the angles are released after being set to 90, gives a different structure than (9), which is the same structure, but where the angles are released from the experimentally found experiment. This suggests there are multiple local minima in the energy landscape and that multiple

conformations are actually possible for this MOF. It is known that MOFs often have very flat energy landscapes and flexible structures; RK26 is likely no exception.

Can we trust simulations to be true to reality?

This raises a further point of discussion: to what degree can we trust the results of the simulation to be true to nature? One important factor to take into consideration is that all calculations done during the simulations assume a temperature of 0K, in order to find the absolute energy minimum. This means that the calculated lowest energy structure might not be the same as it is at ambient temperatures; most likely, the structure is an energy-weighted average over different

conformations. In addition, as mentioned before, the energy landscape of the structure of RK26 is likely very flat, which means that many conformations are possible and a global minimum could be very hard to find. This is evidenced by the fact that the ISIF=4 setting yield different cell angles depending on the starting structure (jobs (9) vs (12)). Finally, one consideration for the accuracy of the simulation is that the ISIF=0 and ISI=4 settings are both fundamentally unphysical by locking certain unit cell parameters to fixed values. Since the ISIF=3 setting could not be used successfully, there is no indication how unlocking the cell volume would influence the resultant structure.

Errors encountered during the simulation

Performing a simulation is not simply pressing a button and waiting for the results; along the way, several errors and crashes were encountered. These came in two variants: the convergence error and the optimisation failure. To understand these, it is important to understand the way a quantum optimisation works. A simulation is performed in several steps. Each step, called an ion step, calculates the new position of each atom in the simulation to create a new structure. It moves the atoms based on the energy and forces, which are calculated in a self-consistent field (SCF) cycle. During an SCF cycle, the positions of the nuclei are fixed and, using density functional theory, the electron density around these nuclei is optimised. In each step of the SCF cycle, the electron density is shifted to lower the energy until it remains constant (self-consistent). For every ion step, about 10-20 SCF steps are performed.

A convergence error happens when during an SCF cycle, the energy does not become self-consistent. This is often the case for MOFs, since their flat energy landscapes often cause the simulation to make large ‘leaps’ across this landscape in order to find the energy minimum quickly. This can cause the simulation to enter an area with an unexpectedly high energy, therefore failing to converge. The convergence error showed up very often; most of the time, telling the simulation to rerun from the last successful ion step was enough to get it running again. If the error showed up more regularly, the simulation could be allowed to continue by setting the EDIFF setting to a lower number. EDIFF determines the break condition for an SCF cycle; when the energy difference between two steps of the cycle is smaller than EDIFF, the cycle stops and proceeds to the next ion step. Reducing EDIFF increases the accuracy of the SCF cycle and therefore avoids inaccurate forces that push the system to unphysical conformations.

Finally, an optimisation failure occurs when the resulting forces from the SCF cycle cause the simulation to undergo unwanted changes during the ion steps. For example, no simulation could be successfully performed with the ISIF=3 setting, where both the cell angles and the cell volume was released. Without fail, these simulations caused the framework to ‘explode’, the unit cell expanding

6 to lengths of over 100 angstrom within a few cycles. A similar

problem arose with simulation (10); the initial structure of the MOF, which was forced to have cell angles of 90 degrees, was so strained that the molecules of the framework fragmented. An optimisation failure can be either numerical (caused by mistakes in the code leading to unphysical values) or physical (where an unrealistic structure is used as input, leading to collapse of the structure). It is possible that the evacuated structure of RK26 is simply not stable; but this does not apply to the structure with filled pores, which we know to exist. Therefore, a numerical error is expected. One explanation for the rapid expansion during the ESIF=3 simulations is the so-called ‘Pulay stress’ that occurs when the simulated cell changes shape: because the k-vectors of the plane wave basis

set of the DFT calculation are still pointing along the vectors of the original cell shape, the calculated forces are inaccurate.

There are several ways to solve issues with optimisation failure. One is to decrease the step size of the ion steps using the POTIM setting. Decreasing the step size means the atoms move less before the next SCF cycle, therefore avoiding the structure exploding. Another is to manually stretch the unit cells by a few angstrom, increasing the distance between atoms and reducing the forces. The former method worked for simulation (10), but neither helped during the ISIF=3 jobs. Another possible solution could be to first perform a classical optimisation (which is much cheaper and faster than a quantum optimisation) to reduce the worst of the strain on the framework, then switching back to quantum to get the final structure. Due to time constraints, this venue was not pursued.

Structural correlations to predicted pXRD patterns

Since the energy landscape of RK26 is flat, large changes in conformation can occur easily under influence of a change in solvent or temperature. In addition changing the solvent in the pores by another might also cause significant changes in the structure and thus, in the observed pXRD pattern. This hypothesis is supported by experimental data: upon soaking in methanol for several days, the measured pXRD pattern of RK26

showed significant changes. Because of this, it is hard to determine experimentally whether a solvent exchange or evacuation experiment was successful. The results of the simulations performed during this project can aid in this. Using Mercury, the pXRD patterns of the finished VASP jobs were examined and compared. The predicted pXRD patterns can be found in the appendix.

All structures have a large peak at =7 that is also seen in many other MOFs. In some structures, this peak is split into two peaks; it is unknown what causes this, but it is likely there is symmetry involved. One trend that can be seen is that the solvent-free structures have the least additional peaks after this first large peak, even less than the completely evacuated structures. The pattern of the clean structures resemble those of the filled MOF, but lacks most of the ‘grass’ present between =15 and =30. These many peaks must then be the result of the solvent

molecules in the pores. Looking at the predicted patterns for the cells with angles of 90°, some large

Figure 3: when cell volume is released (ISIF=3), the simulation caused the framework to 'explode'.

5 15 25 35 45

7 peaks appear in the pattern. It is likely that the increased system of the 90° system cause some smaller peaks to fuse into these larger spikes. As soon as the angles are released in (12), the peaks once again collapse into many smaller peaks.

It is interesting to determine what part of the pXRD patterns is determined solely by the metal atoms of the framework and, by comparison, what is changed by the organic molecules of the framework. Using computational chemistry, this is easily done: after removing all non-terbium atoms from the structure, the pXRD pattern can be predicted by Mercury just as before. Generally, most peaks are still present when the organic molecules are stripped from the metals. Almost every peak increases in intensity relative to that of the first peak at =7 (or, indeed, the peak at =7 decreases in intensity compared to the others). Interestingly, when the process is reversed and the metals are removed, the pXRD pattern looks nearly identical to the structure with metals; only the small ‘grass’ peaks are removed. The spectrum with all atoms included is, with the exception of a few peaks, the sum of the two derivative spectra. It appears that the strongest peaks in a pXRD pattern are caused by both the organic atoms and the metals, while most of the noise is created by the metal atoms. This is the reverse of my expectations.

One problem with the synthesis of RK26 is that judging from the measured pXRD patterns, often structures are formed that might or might not be the same; the measured pattern often differs from sample to sample. Now that simulations have shown that RK26 has a very flexible structure and that the presence or lack of solvents is enough to

radically change the pXRD pattern, some light might finally be shed on why this happens. Several of the predicted patterns more or less resemble deviant experimental data. The same is true for RK26 after solvent exchange. It could be that all along, the right structure was formed, but changes in the contained solvents (which could be caused by a change in solvent concentration, the time a sample has been dry for or even the air humidity) changed the pXRD, causing confusion.

Interestingly, the predicted pXRD of job (11), with a fully evacuated framework and angles locked to 90°, most resembles the pXRD pattern that was found experimentally most often, even more so than that of job (10), which has coordinated solvents. It is difficult to explain why the system at 90° resembles the experimental data the most. The most likely explanation is that because

5 15 25 35 45

Figure 6: One measured pXRD of RK26 'gone wrong'. The pattern resembles the predicted pattern of job (7)

5 25 45

5 25 45 5 25 45

8 the system is so flexible, the pattern (which is measured at room temperature!) comes from an averaged signal from several superimposed conformations, which average out on 90°. To confirm this hypothesis, experiments would need to be performed. A start would be to perform another powder XRD experiment at very low temperatures and seeing if the resulting pattern changes.

Conclusion

Computational chemistry can give valuable insight into the chemistry of solid materials and allows us to study systems that cannot be studied experimentally. After subjecting it to multiple simulations, much insight has been gained into the MOF RK26. The structure is likely very flexible and has a very flat energy landscape. Because of this, multiple conformations are likely at room temperature and large shifts in cell angles can be expected when the solvent in the pores, both free and coordinated, is changed. Using this insight, some questions could be answered that had arisen from experiments. It is likely that perceived ‘mistakes’ in synthesising RK26 had likely been successful, but small changes in the pore occupation caused a large shift in the measured powder XRD pattern, causing confusion. Similarly, a solvent exchange that was thought to have caused framework collapse could actually have been successful as well. With these new insights, the behaviour of this structure can be better predicted, which can cause a boost in the experimental research.

Evaluation and reflection

During the project, I learned several technical skills: I learned to work in a LINUX environment, an operating system I barely had any experience with before, but which is widely used in computational science. During the project, I have used several different programs for the visualisation and editing of structures: Chimera, Vesta and Mercury. Finally, I have gained experience with working with VASP and with performing computational studies in general.

Over the course of the project, I have come to appreciate the essential role of computational models and simulations in chemical research. Doing simulations allows you to study systems and predict their properties in a way that is impossible for experiments. Computational research doesn’t substitute doing experiments, but it does compliment it in order to gain better insight and

understanding in the molecular world. I certainly think that the skills I have accrued over the past two months will come in useful in the future. Because of this, I think this project has helped in making me a better scientists overall.

9

Appendix A: optimised structures

Figure 7-9: the resulting structures of the simulations with ISIF=0 (cell angles locked to experimental values)

Figure 10-12: the resulting structures of the simulations with ISIF=4 (cell angles released).

(1) (9) (8) (7) (3) (2)

10

Figure 13-15: the resulting structures of the simulations with cell angles set to 90° (left, middle) and the following simulation where cell angles were released again (right).

(12) (11)

11

Appendix B: predicted pXRD patterns

Figure 16-18: the predicted pXRD of the simulations with ISIF=0 (cell angles released).

Figure 19-21: the predicted pXRD of the simulations with ISIF=4 (cell angles released).

Figure 22-24: the predicted pXRD of the simulations with cell angles set to 90° (left, middle) and the following simulation where cell angles were released again (right).

5 25 45 (9) RK26cleanISIF4 5 25 45 (8) RK26solventfreeISIF4 5 25 45 (12) RK26cleananglesreleased 5 25 45 (10) RK26anglesto90 5 25 45 (3) RK26cleanISIF0 5 25 45 (2) RK26solventfreeISIF0 5 25 45 (1) RK26ISIF0 5 25 45 (11) RK26cleananglesto90 5 25 45 (7) RK26ISIF4