Quasi-classical modeling of molecular quantum-dot cellular automata multidriver gates

Molecular quantum-dot cellular automata (mQCA) has received considerable attention in nanoscience. Unlike the current-based molecular switches, where the digital data is represented by the on/off states of the switches, in mQCA devices, binary information is encoded in charge configuration within molecular redox centers. The mQCA paradigm allows high device density and ultra-low power consumption. Digital mQCA gates are the building blocks of circuits in this paradigm. Design and analysis of these gates require quantum chemical calculations, which are demanding in computer time and memory. Therefore, developing simple models to probe mQCA gates is of paramount importance. We derive a semi-classical model to study the steady-state output polarization of mQCA multidriver gates, directly from the two-state approximation in electron transfer theory. The accuracy and validity of this model are analyzed using full quantum chemistry calculations. A complete set of logic gates, including inverters and minority voters, are implemented to provide an appropriate test bench in the two-dot mQCA regime. We also briefly discuss how the QCADesigner tool could find its application in simulation of mQCA devices.


Background
Recent advances in molecular electronics on one hand and the limitations of conventional semiconductor devices, on the other hand, have driven a surge of activities towards the realization of molecular devices, circuits, and systems. Achieving the ultimate diminution in size, power consumption, and delay of electronic devices and systems has always been a challenging endeavor of scientists and designers in this field. Due to the prospect of size reduction in electronics offered by molecular-level control of properties, molecular electronics provides means to extend the Moore's law beyond the foreseen limits of small-scale conventional silicon integrated circuits. The small size of molecules allows high device density in the range of 10 11 to 10 12 devices/cm 2 [1]. Besides, the chemical self-assembly capacity in manufacturing molecular devices provides many advantages to conventional semiconductor manufacturing technology, including lower manufacturing cost and uniform device reproducibility. Molecular electronics endeavors to use the nonlinear current-voltage characteristics of individual molecules or molecular assemblies as active devices (diodes, transistors, etc.) in electronic circuits. However, the power consumption of molecular current switches at very high frequencies is still a drawback [2]. The π-σ-π mixed-valence type molecules, which provide double-well potentials for electrons, have been proposed and studied by Aviram towards the synthesis of memory, logic, and amplification [3]. Lent proposed using molecules in representing binary information within the molecular quantum-dot cellular automata (mQCA) paradigm [1,4]. Molecular QCA provides an alternative approach to represent and process data, where binary representation lies in the charge configuration within molecules rather than in the on/ off states of current switches. A cell in the mQCA model consists of a number of molecular quantum dots (or redox centers) and a few electrons. The electrons tend to occupy antipodal sites as a result of Coulomb repulsion. The Columbic interactions cause electrons to tunnel from one redox center to another in a cell, but not between cells. Thus, it is likely that no current flows, since the neighboring cells are coupled by electrostatic field. Figure 1 depicts two-dot and four-dot mQCA cells and how binary "1" or "0" is represented. The first QCA device was implemented and tested using metal dots at near 0 K [5]. Semiconductor implementation of QCA using GaAs/AlGaAs heterostructure materials has been reported in [6,7] as well. Molecules are good containers for keeping electric charge and mQCA cells have a more promising future to work at room temperature [8]. Nonbonding π or d orbitals of a single molecule (or multiple molecules) can function as quantum dots, where the electric charge is localized in each cell. Synthesis of two-dot and four-dot mixedvalence candidate molecules for mQCA has been reported in [9][10][11][12][13][14]. Many of these molecules are mixedvalence type and include transition metals to enable fast electron transfer reactions [15]. Molecular QCA gates are the building blocks of circuits in this paradigm. Calculation of the electronic structure of mQCA gates composed of these molecules is challenging, since the number of basis set functions grows exponentially as the number of molecules and atoms are increased. Besides, many of the ab initio methods fail in describing charge distribution in mixed-valence complexes. Therefore, developing semi-classical models to study mQCA gates is of high importance. Currently QCADesigner [16] utilizes nonlinear and two-state approximations to solve metallic-based QCA circuits. Several QCA circuits including combinational as well as sequential circuits have been studied using QCADesigner. Examples are adders, shift registers, RAM, digital data storage, and simple microprocessors [17][18][19][20][21][22][23][24][25][26]. In this paper, we center on two-dot mQCA and analyze the validity and accuracy of the two-state model approximation for studying multidriver mQCA gates. This study provides an approach to enhance the QCADesigner tool for simulation of mQCA devices in the future.

Two-dot molecular QCA test bench
The majority voter (MV) and the inverter (INV) gates [25] are the fundamental building blocks of any circuit in the four-dot QCA architecture. These gates have been schematically shown in Figure 2a, b. Particularly, the MV gate is referred to as the universal QCA gate, since the AND and OR logical operations can be done by this gate, as evident from the truth table shown in Figure 2c.
Our multidriver minority voter (MinV) gate is composed of m drivers, where m is an odd number, as inputs and one output. Figure 3a schematically illustrates the three-driver MinV model gate in the two-dot mQCA regime. When only one driver (e.g., the d 1 ) is present, the model gate serves as an INV gate (Figure 3b). In the multidriver MinV model gate, all the distances between the centers of the molecules are l, which is equal to the distance between the middle of upper and lower πbonds as shown in Figure 3c. The inputs of the gates are kept fixed, while the output cells switch to their stable states. To this end, two point charges q and 1-q separated by distance l are used to mimic each input as depicted in Figure 3c.
The MinV gate is an alternative to the MV gate in the four-dot QCA, where the output is inverted. Compared to MV and INV gates in the four-dot architecture, which require 16 and 28 quantum-dots correspondingly, the MinV and INV gates require only 8 and 4 quantum-dots in the two-dot mQCA regime. Consequently, these gates provide a small two-dot mQCA test bench, which make high level quantum chemical calculations feasible. The MinV gate can perform NAND and NOR logical operations, as shown in Figure 3d, and provides a functionally complete logic set to implement any logic function in the two-dot mQCA framework. Additionally, it is possible to implement multi-input (or multidriver) MinV gates, which in turn decrease the total number of gates required to implement a logic circuit. It is important to note that since the MinV gate is not a planar gate, circuits implemented in the two-dot mQCA regime are not planar circuits. We highlight that the practical QCA circuits require clocked-control cells and clocking schemes [21,[27][28][29], which are not addressed in this paper.

Two-state model for molecular QCA gates
The charge configuration in a QCA cell is quantified by the so called 'polarization' , and is defined as [30] where q A , q B , q A′ and q B′ are the charges localized at four quantum dots labeled in Figure 1. For two-dot mQCA cells, the polarization is given by the charges q A and q B at the corresponding redox centers in Equation 1. The polarization of a QCA cell varies between −1 and 1, while negative and positive polarizations represent binary "0" and "1", respectively. In two-dot mQCA cells, the normalized dipole moment of the used two-dot molecule is also identical to the polarization, which is given by where μ α denotes the component of the molecular dipole moment that is parallel to the bridge direction, and the origin is in the middle of the bridge. The dipole moment of an mQCA cell can be obtained through full quantum chemical calculations. An important parameter of a QCA device is the Kink energy (E k ), which is the required energy to excite the system from the ground state to the first excited state. To distinguish a bit value from the thermal environment, E k must be greater than k B T [31], where T is the operation temperature in Kelvin, and k B is the Boltzmann's constant. The E k represents the energy cost of cells i and j having opposite polarizations. That is, the electrostatic interaction between all the charges in cells i and j is calculated by [16] where E 0 is the permittivity of free space, and E r is the relative permittivity of the material system. The Kink energy is then given by [16] where E′ i,j and E i,j denote the electrostatic energy of cells i and j having opposite and same polarizations correspondingly. Tougaw and Lent have used a simple Hamiltonian of the extended-Hubbard type to describe the dynamic behavior of four-dot metallic-based QCA nanodevices [32].
Although this Hamiltonian describes the dynamics of the coherent system composed of arrays of four-dot QCA cells elegantly in theory, it is only possible to model the small systems employing this scheme, since the total required number of direct-product basis sets grows exponentially with the number of cells. In other words, an array with N number of four-dot cells and B number of basis sets in each cell requires the total number of direct-product basis sets as [32] By reducing the number of basis sets for each cell and picking up the two orthogonal ones, the Hamiltonian of a four-dot QCA wire can be mapped to Ising model as [32,33] where E k is the kink energy of four-dot cells, γ is the tunneling energy, and σ x and σ z are Pauli spin matrices. In semiconductor and metal-dot QCA, the tunneling barriers of the cells are connected to electrodes, and their heights are controlled externally by voltage sources [33]. The steady-state polarization of any cell, j in a block of cells, can be obtained as a solution to the Hartree-Fock intercellular approximation. This approximation decouples the line of N cells into N single cell subsystems and assumes that the cells are only coupled through expectation values of polarizations. The consequent solution is [33], where P j is the sum of the polarizations of the neighboring four-dot QCA cells. Equation 7 is currently used in the nonlinear and two-state simulation engine of QCA-Designer to solve the metallic-based QCA circuits. It is important to note that mQCA utilizes non-abrupt clocking to reduce the probability of Kink, the property that is not currently present in the QCADesigner as it is based on metallic QCA. In mQCA, the tunneling barriers can be controlled by external electric field [27]. It is demanding to enhance the tool to be able to simulate mQCA circuits. As a primary step towards this end, we present how a similar equation to (7) can be derived directly from the two-state approximation in electron transfer theory [34,35] for two-dot mQCA. We then discuss how these approximations affect the results compared to those obtained from full quantum chemistry calculations. Equation 8 describes the electron transfer (ET) process in a two-dot mQCA cell, where the two redox centers, A and B, are linked through an intervening bridge, I.
The electronic coupling between the redox centers, which is time independent, is an important factor in the ET process. Within the two-state approximation, the Landau-Zener model [36,37] for avoided crossing of energy surfaces may be applied, where the two diabatic states "1" and "0" denoted by ψ a and ψ b , and with energies H aa and H bb interact. The ground state ψ 1 and the first excited state of a QCA cell, ψ 2 , can be related to the diabatic states ψ a and ψ b by a unitary transformation [34] In Equations 9 and 10, ψ 1 and ψ 2 are orthonormal, whereas ψ a and ψ b are not orthogonal in general. The correspondence between diabatic and adiabatic two-state models arises from the secular determinant (S ab = <ψ a |ψ b > is neglected) [38] where E is the adiabatic energy eigenvalue. The η in Equations 9 and 10 satisfies [38] tan2η The energy difference between the two diabatic states in the output cell of the MinV gate can be approximated by calculating the difference between the electrostatic energies of the gate for the two output configurations, where the unit charge is localized at sites A and B correspondingly. Using Equation 3, for the "1" and "0" output states (Figure 3a), we obtain Inserting Equation 13 into Equation 12 we have The Kink energy of two-dot cells can be calculated from Equation 3 and 4 for two neighboring cells as and for each driver, the polarization is defined using Equation 1 Thus, Equation 14 can be rewritten in terms of the Kink energy and the input polarizations And finally, using Equation 1 and the transformation coefficients in Equations 9 and 10, the output polarization of the MinV gate is obtained Inserting Equation 17 into Equation 18, we can find the polarization of the output cell of the MinV gate in terms of the polarizations of the inputs straightforwardly as Equation 19 in two-dot mQCA is analogous to Equation 7 in four-dot metal-based QCA, where the tunneling energy γ appears as electronic coupling of redox centers (H ab ) in Equation 19. They also imply The additivity relation in Equation 20 originates from the additivity of electrostatic potential energy in Equation 13 for diabatic states.
Multidriver MinV gates help reduce the number of needed gates for implementation of a logic circuit. Similarly, for an m-input MinV gate, we obtain We refer to Equaiton 21 as the two-state model (TSM) for mQCA gates along this paper. The E k and the H ab are the only parameters of the TSM. Once the geometrical parameter l is determined, experimentally or from theoretical calculations, the Kink energy can be calculated using Equation 15. The electronic coupling matrix element, H ab , can be calculated using various quantum chemistry techniques [34,[38][39][40][41] or obtained via spectroscopic experiments, including absorption [42,43], EPR [44], and ultraviolet photoelectron spectroscopy [45]. As we will present, the parameter μ = E k /2H ab plays an important role in the accuracy of the TSM. It is also the slope of the switching response function at the origin i.e., Results and discussion The chemistry of mixed-valence complexes has received considerable attention recently in mQCA device implementation, where the intramolecular electron transfer and charge localization at redox sites are the important key factors. Mixed-valence compounds contain more than one redox site in the same molecule or molecular unit. Simple model molecules for two-dot cells are the πσ-π mixed-valence types, which were proposed by Aviram and studied later by Hush [3,46]. In the Aviram's model molecule (1, 4-diallyl butane cation), the two πallyl groups form two redox centers and are connected by a σ-butyl bridge. One of the allyl groups is a neutral radical, while the other one is anionic (or cationic). The possibility of charge localization in some mixed-valence mQCA candidate molecules has been examined theoretically as well as experimentally [15,47,48]. Advances in quantum chemistry in the past half century provide reliable methods to explore the electronic structure of molecules; however, many of the ab initio techniques fail in describing charge distribution in mixed-valence complexes. The unrestricted Hartree-Fock method overestimates the charge localization due to the neglect of electron correlation effects [49]. In the density functional theory (DFT) method, the exchange potential defined in hybrid functional leads to underestimation of charge localization [47,49]. The complete active space selfconsistent field (CASSCF) method [50,51] is believed to be the most reliable for describing charge distribution in mixed-valence complexes [50]. However, the multideterminant CASSCF calculations scale with the system size, which makes this method highly demanding in computer time and memory. The number of Slater determinants has factorial dependence on both the number of active electrons and particularly on the number of active orbitals generating many-electron configurations (full configuration interaction (CI) within the active space). This is much more significant than any dependence on the number of one-electron basis functions. The number of Slater determinants in a full CI calculation is given by: where M is the number of active orbitals, N α and N β are the numbers of active electrons with αand β-spins, respectively, and the quantities in parentheses are binomial coefficients: We present full quantum chemistry calculations of the steady-state output polarization of the universal MinV model gate serving as INV and three-input MinV gates. The results based on full quantum chemical calculations are compared to the results obtained from the TSM. The π-σ-π mixed-valence type molecules, descended from Aviram's original idea are analyzed. These molecules include 1, 6-heptadiene, 1, 8-nonadiene, and 1, 4-diallyl butane radical cations, which will be referred to as molecule 1, molecule 2, and molecule 3 in this paper, respectively ( Figure 4). We optimized the geometry of these monocations using the DFT/B3LYP method. The dot-dot distance, l, in these molecules is between 0.5 to 0.8 Å. Bistability and electron localizability of these molecules have been studied in [3,46,49].
Koopmans' theorem [52] has found extensive application in calculation of the ET matrix element, H ab for symmetric molecules. Under the two-state approximation, H ab is related to adiabatic energies of the ground and first excited state (E 1 and E 2 ) as [34] When no driver is present or the sum of the input drivers is zero, P o = 0; and from Equation 18 it is clear that Cos2η = 0, thus According to the one electron Koopmans' theorem, the ionization potential of the highest occupied molecular orbital (HOMO) and HOMO-1 can be expressed in terms of the molecular orbital (MO) energies, i.e., [38][39][40][41] I HOMO ¼ ÀE HOMO ð27Þ Inserting Equations 27 and 28 into Equation 26, the electronic coupling element is obtained in terms of the MO energies as [38][39][40][41] The state-averaged CASSCF (SA/CASSCF) method [53] can be used to calculate the electronic coupling element Figure 4 Geometry of the molecules we used in our calculations. 1, 6-heptadiene, 1, 8-nonadiene, and 1, 4-diallyl butane radical monocations from left to right. In the first two molecules, coordinates of the three highlighted atoms have been set to xy plane, while the coordinates of the central carbon have been set to origin. For 1, 4-diallyl butane, the origin has been set in the middle of the bond between the two central carbon atoms. of asymmetric molecules. One can obtain H ab by calculating the energies of the ground and first excited states and use them in Equation 26 within the two-state approximation. We have calculated the electronic coupling elements using both methods. The calculations for molecule 1 and 2 were based on SA/CASSCF (3,4), and the calculations for molecule 3 were based on SA/CASSCF (5,6). In 1, 4diallyl butane cation, the allyl π-bonds are aromatic, and the active space is extended to five electrons in six orbitals. The calculated ET matrix elements have been compiled in Table 1. The Kink energies of the molecules have also been calculated using Equation 15. Table 1 includes all required parameters of the TSM. All calculations reported here were performed in Gaussian 09 [54], using 6-31 G(d) basis set. Various basis sets have been extensively tested to examine the basis set dependency of the results. Application of larger basis sets did not significantly influence the energy difference between the ground state and the first excited state. The results from the two methods are in good agreement. The ET matrix element, H ab decays exponentially with dot-dot distance, l [55]. The dot-dot distance for molecule 3 is less than that of molecule 2; however, the H ab has been remarkably decreased. This is due to the symmetry of this molecule and the aromatic bonds of the radical allyls. The geometrical parameter, l, and the type of head groups, play an important role in determining the ET matrix element. To obtain a more accurate electronic coupling element, the overlap integral, S ab , should be taken into account as described in the work of Farazdel et al [56]. Aviram [3] has obtained a negligible overlap integral for molecule 3 with dot-dot distance of 7 Å. In mQCA, the electron transfer drama should have a little effect on the geometric parameters [4]. Consequently, candidate molecules should possess fast electron transfer reactions, and the relaxation of nuclear degrees of freedom can be ignored. Table 1 also lists the changes in the head groups' π-bonds (Δζ) as a consequence of ET reactions. It is seen that ET reactions in molecule 3 should be faster compared to the other two molecules.

INV gates
The INV gate in two-dot mQCA is the nucleus of all other gates. Once its operation and switching properties are clearly understood, the properties of more intricate structures such as multidriver MinV gates can be derived from extrapolating the results obtained from the inverters, based on the additivity relation (Equation 20). The analysis of inverters can be extended to explain the behavior of more complex gates, which in turn form the building blocks for modules such as adders, multipliers, and processors. Table 2 compares the output polarizations of the INV gates, obtained from full quantum chemistry calculations and the TSM. The normalized dipole moments (Equation 2) of the monocations adjacent to fixed inputs (point charges as fixed drivers) have been calculated based on SA/ CASSCF(3,4) for molecule 1 and 2, and SA/CASSCF (5,6) for molecule 3. The root mean square errors (RMSE) of the results obtained from the two methods have been calculated. The RMSE decreases with the increase of the μ parameter (or decrease of the ET matrix element, H ab ), determining the degree of agreement between the results. The saturation polarization of the output is also dependant on the μ parameter. For INV gates, it is obtained by setting the sum of the input drivers to one in Equation 21 as Equation 30 shows that the saturation polarization of the output increases with the increase of μ. This is also evident from the results in Table 2.  When there is a single driver, the change in the sign of the input will result in a change in the sign of output polarization. The negative inputs have been omitted to reduce the size of the table. P o and P o * denote the calculated output polarizations based on the SA/CASSCF and TSM methods respectively. RMSE * represents the root mean square errors of P o and P o * .

Two-driver devices
In the model MinV gate, the number of input drivers (m) should be odd. No logic operation is performed when m is even. However, neglecting the logic, the twodriver device is an appropriate small model system for studying the additivity relation, and how the accuracy of the TSM is influenced by the number of drivers. Here, the MinV gate is probed when only the two input drivers d 1 and d 3 are present (Figure 3a). We have calculated the normalized dipole moments of the gates' outputs based on the SA/CASSCF calculations. The output polarizations have also been calculated by the TSM. The results obtained from the two methods are compiled in Table 3, which are in good agreement. The conclusions from analysis of the INV gates can be extrapolated to two-driver devices as well. Compared to INV gates, the increase in the RMSE of the two-driver devices, composed of molecule 1 or 2, is mainly attributed to the asymmetric head groups. In other words, the effect of d 1 on the head groups is different from that of d 3 , where P d1 = P d3 . In molecule 3, the allyl head groups are symmetric, and the TSM error mainly arises from the classical approximation of the intercellular interactions. It is important to note that the output polarization of the two-driver devices can be calculated by employing the additivity relation on the output polarizations of the INV gates. The additivity relation has been validated for the SA/CASSCF method as well. Through full quantum chemistry calculations, the output polarizations of the two-driver devices were obtained. We also used the results of the INV gates, P(P d1 + P d2 ) from Table 2, to examine the additivity relation for SA/CASSCF method. As expected, RMSE is highly dependent on the symmetry of the head groups. Unlike molecule 1 and 2, for the case of molecule 3, each driver has an exactly same effect on the head allyl groups, which leads to smaller RMSE. We also highlight that employing additivity relation decreases the computational cost of SA/CASSCF calculations. Table 2 and Table 3 also show how the accuracy of the TSM is affected by the number of input drivers. It is seen that RMSE is approximately doubled when the number of input drivers is scaled up by a factor of two.

Three-input MinV gates
The multidriver MinV gate is a universal gate in the two-dot mQCA. The output polarizations of these gates with three fixed input drivers are shown in Table 4. This table also shows that the output polarizations of the MinV gates can be obtained from extrapolating the INV output polarizations using the additivity relation. For the model molecules, considering the spatial location of the d 2, the effect of d 2 on the head groups is different from the same effect of d 1 and d 3 , while in the TSM, drivers with equal polarizations have same effects on the head groups and are treated the same. Quantum chemical calculations show that despite the equal sum of the input polarizations, the output polarizations are not equal particularly when the sum of the drivers is zero. Table 4 also displays that the SA/CASSCF method returns different output polarizations, while the sum of the input polarizations is zero. This is the main reason of the decrease in the accuracy of the results obtained from the TSM for MinV gates. Ignoring these points by avoiding the null state logic occurrence, the two-state approximation results are fairly in good agreement with the quantum chemical calculations. Table 2 and Table 4 show how the accuracy of the two-state model is decreased with the number of drivers. It is seen that RMSE is tripled when the number of input drivers is scaled up by a factor of three.  The degree of agreement between the TSM and quantum chemical calculations is highly dependent on the μ parameter and the symmetry of the head groups. Additionally, application of the additivity relation for CASSCF method can in turn reduce the computational cost. It is important to note that we did not address questions of surface attachment, input/output, clocked control, layout, and patterning, which are the requirements of a practical QCA system. Moreover, we did not consider the relaxation of nuclear degrees of freedom associated with electron transfer. It is presented that for mQCA, the electron localization and Coulombic interactions play the key roles, and nuclear positions can be considered frozen (nuclear relaxation even assists charge localization) [4]. Although we limited our focus on the two-dot mQCA, it merits highlighting that the model can also be used for four-dot cells, since they can be considered as double two-dot cells. Our focus was on the mQCA gates as building blocks of circuits. The two-state model may be applied to simulate mQCA circuits as well, as it is currently used iteratively for simulation of metallic QCA circuits in the QCADesigner. However, to determine the additive error resulting from exploiting the two-state model for solving mQCA circuits, further quantum chemical calculations on the mQCA clocked circuits composed of several molecules are required, which are extremely challenging at the time, and have not been addressed in this paper.

Competing interests
The authors declare that they have no competing interests.
Authors' contributions ER developed the theories and carried out the quantum chemical calculations. SMN supervised the project. All authors read and approved the final manuscript. P o and P o * denote the calculated output polarizations based on the SA/CASSCF and TSM methods, respectively. RMSE * has been calculated based on P o (P d1 , P d2 , P d3 ) and P o * (P d1 + P d2 + P d3 ). RMSE ** has been recalculated when points with P o * = 0 have been omitted.