Computational analysis to study the effect of infusion of Tetracyanoquinodimethane in zinc based metal-organic framework

The present study reports the computational analysis of band structures and density of state (DOS) profiles of a Tetracyanoquinodimethane (TCNQ) infused IRMOF-1 (Iso Reticular Metal-Organic Framework-1). IRMOF-1 is a zinc-based MOF and has been reported useful in literature for a variety of applications, including gas storage and sensing owing to its highly porous structure. We have adopted a SIESTA based computational investigation approach to investigate the tuning of the band structure and DOS of IRMOF-1 via the infusion of a redox active species, i.e., TCNQ. The results have highlighted that the doping of TCNQ in IRMOF-1 is useful for the realization of a novel material with an estimated band gap of 4.35 eV. The necessary computational steps involved various optimization steps, such as the optimization of mesh-cutoff, lattice constant, lattice volume, and conjugate gradient. These computation studies have thus established that an approach of doping IRMOF-1 with TCNQ can result in the development of a suitably functional porous material that can be further exploited as an ultrawide‐bandgap semiconductor material or for the electrochemical sensing of different analytes, such as gases.


Introduction
Metal-organic frameworks (MOFs) are the nanoporous materials composed of metal clusters connected with organic linkers. MOFs have attracted significant interest in the last decade for a variety of applications [1]. MOFs offer wide ranges of pore sizes, permanent porosities, large surface areas, and chemical/thermal stabilities. The structural tunability of MOFs allows their modification with other suitable materials (e.g., via doping) so that the resulting composite structures become useful in diverse applications [2]. With respect to the range of applications, MOFs have found profound usage in various chemical and biological areas, including gas storage and separation [3], molecular/environmental sensing, catalysis [4], drug storage, and delivery [5].
For the last two decades, the first-principles simulation of condensed matter systems has expanded spectacularly, from physics and chemistry into life, earth, nanoscale and materials sciences [6]. This success has been driven by both the steady growth of hardware computing power and the software development. The importance of hardware resource development is directly realized for the analysis of many deeply rooted algorithms within the standard methods as the researchers need to process computation with increased number of atoms in the simulation box (scaling as N 3 ). Though cube-scaling algorithms are still used in most of the DFT (density-functional theory) calculations (up to about 1000 atoms on the best supercomputers), it is obvious that linear-scaling algorithms will be of advantage in the long run [6][7][8].
The computational analysis of various parameters of MOFs has been mainly based on first-principles DFT within the Generalized Gradient Approximation (GGA) of Perdew-Burke-Ernzerhof (PBE) functional. DFT methods have been in use for several years to study the electronic structures of many-body systems. The DFT roots can be traced as far back as the mid-1960s, in Hohenberg-Kohn theorems (H-K). In DFT, 'Functional' Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. simply refers to the function of a function. The exchange-correlation functional is set to Generalized Gradient Approximation (GGA, PBE version) and the basis size is set to Double Zeta with Polarization (DZP). Conjugate Gradient (CG) minimization algorithm is used to perform structural relaxation. These methods take in account the parameters like numbers and types of atoms, molecules, and condensed phases. In principle, main objective of the above studies has been to assess the ground state of systems under investigation.
The practical quantum mechanical simulations of material, which conventionally takes in account the electronic degree of freedom, are limited to only about 1000 atoms at one time [9]. In contrast, the classical simulations employing empirical interatomic potentials may involve the computation to over 10 9 atoms. Much of this 10 6 -factor difference is due to the existence of the well-developed order-n algorithms for the classical methods, in which the computer time and memory scale linearly with the n-number of atoms of the simulated system [9,10]. The major task of any computer simulation in the area of material science is to predict the behavior of materials under varying conditions.
The computer aided simulation studies provide an important tool to the researchers to develop and personalize the new materials vis-a-vis an intended application. Particularly in context to the advancement in nanotechnology, the computer simulation studies on nanomaterials have assumed a great significance as the nano dimensional materials show much different properties than their bulk counterparts. The popularly used computational software for molecular simulation of materials and the visualization of the results thereof include Gaussian and Gaussian View, CP2K, VASP, and SIESTA. Gaussian can simulate the materials' behaviour based on the predictions of energies, molecular structures, vibrational frequencies, and molecular properties [11]. CP2K is a free-available software package, written in Fortran, and distributed under the GPL license [12]. It can perform atomistic simulations of solid state, liquid, molecular, periodic, material, crystal, and biological systems [13]. VASP (Vienna Ab-initio Simulation Package) is a commercially available computer program for performing atomic scale materials modelling, e.g., electronic structure calculations and quantum-mechanical molecular dynamics [14]. VASP can run on super-scalar processors, vector computers, and parallel computers. SIESTA (Spanish Initiative for Electronic Simulations with Thousands of Atoms) consists of method(s) libraries and functions toward their computer program implementation [15,16]. SIESTA is used to perform electronic structure calculations and ab initio molecular dynamics simulations of molecules and solids [17]. This software is written in Fortran with provision of dynamically allocated memory. It may be compiled for serial or parallel execution (under MPI). SIESTA is free for academic use and has the following computational capabilities: Total and partial energies, Atomic forces, Stress tensor, Electric dipole moment, Electron density, Geometry relaxation, Band structure, Atomic, orbital and bond populations [17].
In the present computation study, we have employed SIESTA based computational approach for the assessment of band structure of a zinc-based MOF, i.e., IRMOF-1, before and after its doping with an external organic molecule 'Tetracyanoquinodimethane (TCNQ)'. IRMOF-1 has been selected for the study as this was one of the first MOFs synthesized and is widely regarded as a prototypical MOF with simple topology. IRMOF-1 bears a low framework density and significant amount of void spaces. The advantages of facile synthesis and availability of void space make it an attractive option with respect to many technological applications, including gas storage and gas sensing. On the other hand, TCNQ is well-known to improve the electronic properties of certain other MOFs, e.g., Cu-BTC [8]. We envisage that the TCNQ doped IRMOF-1 can open up an exciting opportunity to expand the applications of this new functional material as a potential ultrawide-bandgap semiconductor candidate or in the electrochemical sensing of various analytes. As such, this article is the first computational report to establish the effect of infusion of TCNQ on the important electronic properties of IRMOF-1, e.g., band structure, density of states, and band gap.

Computational details
Identification or evaluation of the correct ground state of electron density that can achieve an energy minimization is a challenge in the quantum mechanical modeling. The present study has employed SIESTA (Spanish Initiative for Electronic Simulations with Thousands of Atoms) code consisting the framework of Kohn-Sham DFT [7]. The simulations were performed with resources from the High Performance Computing (HPC) facility of the CSIR Fourth Paradigm Institute, Bangalore, India. The SIESTA (version 4.1) was compiled with Intel parallel studio 2017 to run the application in parallel mode. The computation has been applied on a total of 800 number of processes (Number of Nodes=50; Number of CPUs=16; MPI task per node=16), involved in different simulation steps.
A simulation process has four basic attributes: input, simulation block, post-processing, and output. All the SIESTA algorithms are located in the simulation block diagram that process parameters like geometry optimization, molecular dynamics (MD), etc, as depicted in figure 1. For efficient computation with SIESTA, input must be in the correct format and therefore Flexible Data Format (FDF, with extension * .fdf) was used as the main input file.
The present research studies the computational parameters on IRMOF-1 before and after the infusion of TCNQ (Tetracyanoquinodimethane) in its structure. IRMOF-1 ( figure 2) is a Zinc-based MOF, formed by the reaction of Zn 4 O(-CO 2 ) 6 (metal center) with benzene-1,4-dicarboxylic acid (organic linker). It can be easily synthesized and possesses some unique material properties, such as large surface to area ratio, accessible pores, and thermal stability. IRMOF-1 has been reported useful in many applications, including storage of gases, molecular separations, and sensors [18]. Our present computational study attempts to evaluate the possibility of using it as a potential material for future sensor applications. IRMOF-1 is known for its potential to capture many analytes, including gases like CO 2 .
The structure of IRMOF-1 (figure 2), has been obtained from the COD (Crystallography Open Database) as a CIF file. Jmol was used to visualize the CIF file. Jmol console commands were used to export the Cartesian coordinates corresponding to the atoms in the desired .xyz format. Python code was used to translate the Cartesian coordinates to the SIESTA input file format, i.e., .fdf (Flexible File Format). The .fdf file format is the standard SIESTA input file format which can be used to provide the simulation input in any order, while also being able to add other parameters during computations. SIESTA also requires a pseudopotential file for each species to perform the computation. These Pseudopotential files for different species were imported from NNINC (http://nninc.cnf.cornell.edu).

Results and discussion
Computational studies with pristine IRMOF-1 During simulations, the total number of atoms in a molecule needs to be specified. The number of atom species contained in the molecule is also given as an input. For instance, IRMOF-1 has four number of atoms species, i.e., Zn, H, O, and C and the total number of atoms is 424 in the unit cell. SIESTA program functions by matching the values of these parameters with the number of positions it finds in the supplied spatial coordinates. The coordinates are specified in Armstrong unit. Simulations were performed using a single unit cell of IRMOF-1 (a=b=c=25.83200; and α=β=γ=90.0000) with 424 numbers of atoms.
During simulation, selection of suitable mesh of the electron density is an important parameter. The mesh size is also called the mesh cut-off energy. By performing the calculations with different mesh sizes, we computed the single point energy to find out the optimal cut-off energy. As a first step, mesh-cutoff optimization was studied at a periodic interval of 50 starting from 100 to 800 Ry mesh cutoff. Different mesh-cutoff values were extracted from the SIESTA output file during each step of calculation of total ground state energy. The mesh size was tweaked until the computed energy changed only by small amounts. The final total energy [eV] was then plotted against mesh-cutoff [Ry] by an open source software 'Magic-plot' to evaluate the energy profiles of the  After optimization of mesh-cutoff, the lattice constant was optimized within the range of 25.40 Å to 26.75 Å at an interval of 0.05 Å. At all the computed points, a mesh-cutoff of 650 Ry (as optimized and extracted from total energy SIESTA output file) was given as an input parameter. The plot of total energy [eV] versus Lattice Constant [Å] is shown in figure 3(B). The output of this study was used to compute the total volume by the Brich Murnaghan equation [19]. The necessary computations are depicted in figure 4(A). The results have revealed a value of lattice constant as 25.8508 Å which is in a close agreement with previously reported experimental values (25.67−25.89 Å) [20,21].
The final structure optimization of IRMOF-1 has been performed by providing the input parameters of optimized value of mesh-cutoff, Kgrid Monkhorst-Pack, and lattice constant. For computations, the Conjugate Gradient (CG) steps were operated from 1 to 50. The structure could be optimized before the 50 th iteration of the DFT calculation and the related results are shown in figure 4(B). The DFT calculations for band structure and density of states (DOS) have been performed on the optimized/relaxed IRMOF-1 structure. For the assessment of band structure, we have used the seek path [22] tool, which allowed us to find the k-path in first Brillouin zone path on the optimized IRMOF-1 structure. The computation has revealed a Fermi energy value of 4.28 eV. The band structure and DOS are shown in the figures 5(A) and (B). The projected density of states (PDOS) and the total density of states were also obtained. Fermi energy in these graphs is indicated by a horizontal line (in case of band structure based analysis ( figure 5(A))) or a vertical line (in case of PDOS based analysis ( figure 5(B))). As can be seen in figure 5, position of Fermi level toward the valence band indicates a greater number of holes in comparison to electrons, indicating p-type characteristics of the material. The band structure in figure 5 also shows that Fermi level has been shifted to around −4, i.e., toward the valence band. Hence, the material should be exhibiting p-type semiconducting characteristics. In conduction band, major DOS starts around 0 eV while few orbitals before that are occupied. The band gap is estimated to be around 8 eV. It may be outlined here that DFT estimations generally underestimate the band gap.  Computational modelling after diffusion of TCNQ in IRMOF-1 Tetracyanoquinodimethane (TCNQ) is an organic compound possessing an electron acceptor property. TCNQ is useful to forms charge-transfer salts that can offer the development of new materials with high electrical conductivity. Such TCNQ infused materials have been recommended with potential applications in many applications areas, e.g., sensors and organic electronics [23]. Realizing the importance of IRMOF-1 as a novel porous molecule that can offer the development of many types of sensors, including those for gases, we have carried out the computational analysis of the effect of diffusion of TCNQ in IRMOF-1 on band structure and DOS. TCNQ (C 12 H 4 N 4 ) structure has been introduced inside the optimized structure of IRMOF-1. The optimized molecular structures of both IRMOF-1 and TCNQ infused IRMOF-1 are shown in figure 6.
The above structure was optimized by following the steps similar to those adopted in the case of pristine IRMOF-1. The complete structure optimization of TCNQ@IRMOF-1 has been performed after evaluating all the required optimized value of mesh-cutoff, kgrid Monkhorst-Pack, and lattice constant. The structure was optimized in different CG steps from 1 to 105. The progress of steps during the structure optimization within the above interval is shown in figure 8(B).
Finally, the band structure and DOS for the TCNQ@IRMOF-1 has been computed on the optimized structure and Fermi energy value has been found as 5.978 eV. During the process, the different k-path in the first  Brillouin zone were identified using seek-path, as required for SIESTA input. The band structure and DOS are shown in figures 9(A) and (B). In comparison to the band structure of IRMOF-1, a shift in the Fermi level towards the conduction band after the infusion of TCNQ indicates a clear enhancement of conductivity. This is also reflected by a visible shift in the Fermi level (depicted in PDOS diagram). The contribution of N (Nitrogen) atoms in valence as well as in conduction band (CB) has resulted in the narrowing of the band gap. The narrowing of the band gap can be attributed to anti-crossing interactions amongst the localized N states and the extended CB states of IRMOF-1. As the CB of IRMOF-1 actually splits into sub bands, the introduction of N atoms results in the shifting of the lower sub band edge. Having high electronegativity and small size, N atoms   tend to replace heavier atoms of the IRMOF-1 structure from the CB, and, consequently, cause perturbation of crystal lattice potential as evident from the visible Fermi level shift [25]. Still similar to IRMOF-1, TCNQ@IRMOF-1 exhibited p-type nature.
A summary of the present computational study is presented in figure 10. It is clear that the infusion of TCNQ within IRMOF-1 is helpful to tune the bandgap of the material. As such, this bandgap value (4.35 eV) of the TCNQ@IRMOF-1 suggests the potential application of this new material as an ultrawide-bandgap semiconductor material. TCNQ@IRMOF-1 can also be explored for the development of electrochemical sensor for analytes like gases and small molecules.

Conclusions
The SIESTA based computational analysis has been performed to identify the band structures and projected density of states of IRMOF-1, before and after the infusion of TCNQ. The results of both the band structures and projected density of states are in agreement with each other thereby highlighting a robust selection and optimization of all the important parameters during the present computational approach. These studies have revealed that the band gap energy of IRMOF-1 can be tuned by the infusion of TCNQ. As such, the pristine IRMOF-1 has been found to have a band gap of 9.64 eV, which was modified to 4.35 eV after the loading of TCNQ in the structures. This tunable behavior of IRMOF-1 could be of practical significance in the development of amperometric or similar types of electrochemical sensors for a variety of analytes. Thus, the present study should be open up the opportunities toward establishing the usefulness of TCNQ doped IRMOF-1 as an electrochemically active platforms for related applications. Being porous in nature, such IRMOF-1 based sensing platforms can be targeted for a variety of analytes, e.g., gases, anions, volatile organics, etc.