Quantum hardware calculations of periodic systems with partition-measurement symmetry verification: simplified models of hydrogen chain and iron crystals

Running quantum algorithms on real hardware is essential for understanding their strengths and limitations, especially in the noisy intermediate scale quantum (NISQ) era. Herein we focus on the practical aspect of quantum computational calculations of solid-state crystalline materials based on theory developed in our group by using real quantum hardware with a novel noise mitigation technique referred to as partition-measurement symmetry verification, which performs post-selection of shot counts based on $Z_{2}$ and $U_{1}$ symmetry verification. We select two periodic systems with different level of complexity for these calculations. One of them is the distorted hydrogen chain as an example of very simple systems, and the other one is iron crystal in the BCC and FCC phases as it is considered to be inaccessible by using classical computational wavefunction methods. The ground state energies are evaluated based on the translational quantum subspace expansion (TransQSE) method for the hydrogen chain, and periodic boundary condition adapted VQE for our iron models. By applying these techniques for the simplest 2 qubit iron model systems, the correlation energies obtained by the hardware calculations agree with those of the state-vector simulations within $\sim$5 kJ/mol. Although the quantum computational resources used for those experiments are still limited, the techniques applied to obtain our simplified models will be applicable in essentially the same manner to more complicated cases as quantum hardware matures.


I. INTRODUCTION
"Quantum computational chemistry" is expected to be a promising alternative to its conventional (classical) counterpart due to the efficient way in which quantum computers handle large Hilbert spaces.[1,2] Development in this field has been significantly accelerated by the implementation of cloud-accessible quantum computers from several providers.Quantum hardware is still in its early stages, with devices displaying modest numbers of noisy qubits (typically 50-100) that prevent the application of error correction schemes.These near-term devices are frequently referred to as Noisy Intermediate Scale Quantum (NISQ) [3] devices.
Variational algorithms are believed to be the most suitable techniques for NISQ devices by constructing a hybrid quantum-classical setup, in which a relatively shallow parameterized quantum circuit performs heavy tasks such as encoding correlated molecular wavefunctions to calculate the expectation value of the energy, while the classical computer collects the data from the quantum computer to optimize the parameters within the variational loop.The Variational Quantum Eigensolver (VQE) algorithm [4] is one of the most frequently used variational algorithms to run quantum chemistry simulations on NISQ devices.
In the VQE algorithm, we consider a type of cost function E(θ) depending on a set of parameters θ describing an associated trial wavefunction |Ψ(θ) , which is ex-pressed as where Ĥ is a Hamiltonian describing the system.The trial wavefunction |Ψ(θ) is actually implemented into the quantum computer in the following form where |Ψ 0 corresponds to an initial state that should be easily prepared.The unitary operator Û (θ) is implemented as a parameterized quantum circuit, or ansatz.While molecular ground state calculations have been focused on as the target of early-stage quantum simulations [2,5], a number of solid-state quantum calculations have been performed by employing Periodic Boundary Conditions (PBC) [6][7][8][9][10][11][12][13][14].These simulations, however, have been performed either on emulators run on classical machines, or on simplified tight-binding models on quantum hardware.Here, we consider the electronic structure hamiltonian explicitly.Solid state calculations are, in general, more expensive than molecular cases, because the Hamiltonian describing the system depends not only on the spin-orbitals but also on the reciprocal space points (k-points).[15] Most of the efforts devoted to classical computational PBC calculations have been made based on Density Functional Theory (DFT) [16,17].For example, first principles evaluation of the properties of solid-state iron materials has been addressed with DFT methods in many works [18][19][20][21][22][23].Such metallic solid state calculations are known to be sensitive to the choice of DFT exchange-correlation functional without any way to systematically improve the results, in contrast with wavefunction methods where a clear methodology may be followed to obtain increasingly accurate simulations.Using again the iron example, LDA cannot identify the most stable phase of this metal and more complex functionals (GGA in this case) are needed to address this problem [24].Despite the computational difficulty in scaling, wavefunction methods have also been reported [25,[25][26][27][28] to characterize the effect of electronic correlation, which is known to be important in, for example, the interpretation of the behavior of Mott insulators [29].Quantum chemistry calculations on quantum hardware can potentially extend the applicability of wavefunction methods because of the capability of addressing an exponentially growing computational basis with a relatively small number of qubits.
Even at this early stage, it is essential to perform real hardware experiments for quantum computational methods to understand the limitations to be overcome, because actual quantum and/or classical computational costs involved in variational algorithms are not only bound by theoretical computational complexity results [30], but they may also be influenced by the characteristics of the hardware used.Such information should be continuously updated to be utilized for helping hardware development, so that we can achieve quantum advantage as soon as possible.Our motivation to write this paper is framed in this context, where we aim to perform benchmark calculations for representative models in order to have a glimpse of the current state of quantum computing techniques for the study of periodic systems.In particular, we set two targets of PBC calculations as follows: i) distorted hydrogen chain as one of the simplest periodic systems, and ii) iron crystals as a representative complex system that will not be easily accessible by classical computational wavefunction methods.
Quantum hardware development is in its early stages and is probably not mature enough to find quantum advantage in the field of quantum chemistry at this moment.Here instead we demonstrate techniques to reduce the quantum computational resource requirement by providing a series of (approximate) techniques to allow us the execution of quantum hardware calculations with currently available NISQ devices and the study of their performance and challenges in the simulation of periodic systems.The reduced capabilities of current devices force us to use a series of drastic approximations on our models.In order to reduce qubit number requirements and circuit depth, we are restricted to small basis sets, which capture only a fraction of the correlation energy for our systems, and small k-point meshes, which are insufficient to reach the thermodynamic limit.However, these approximations can be rolled back as the hardware improves in the future to scale up the complexity of the target system.
We demonstrate that noise mitigation techniques are essential for obtaining accurate energies from NISQ devices.In particular, we show that State Preparation And Measurement (SPAM) [31] noise mitigation improves the accuracy of the energy obtained from quantum hard-ware, compared to the exact value obtained by statevector simulations on a classical computer.In addition, we introduce a novel noise mitigation technique called Partition-Measurement Symmetry Verification (PMSV) which significantly improves our results.As we will demonstrate later, the combination of these noise mitigation techniques gives us energy values within ∼5 kJ/mol of agreement with theoretical values in the PBC-adapted VQE calculations for our iron models.
We focus on the challenges of calculating total energies in this particular paper, which are important for a wide range of problems including phase stability or defect formation energies.Quantum algorithms for many-body band structure calculations and derived properties will be studied in future work.
This paper is organized as follows.In Sec.II, we provide the computational method of the simulations.In Sec.III, we provide our algorithm benchmarking on simple hydrogen lattices, followed by that of iron crystals starting with classical calculations to determine the best way to simplify the system for calculations using quantum algorithms.Section IV concludes this paper with a summary of our findings and an outlook for the future.

II. METHODS
In this section, we describe the methodology used to perform our experiments on quantum hardware.In order to run these experiments, we follow a pragmatic line and, for each system considered, we choose a quantum algorithm that allows us to create the simplest quantum circuit compatible with the particular features of the system considered and thus ensures that it will be reliably executed on the quantum computer.Following this criterion, Translational Quantum Subspace Expansion (TransQSE) is applied for the hydrogen chain, whereas VQE with unitary coupled cluster singles and doubles with PBC (UCCSD-PBC) ansatz is used for iron crystals.The main gain is in using conservation of crystal momentum to simplify the ansatz, by discarding redundant excitations.This reduction suppress the scaling of the number of parameters in ansatz from O(N 4 L 4 ) to O(N 4 L 3 ), where N and L are the number of orbitals in a primitive cell and the number of k points, respectively.For further details about these quantum computational theories, see Ref. [13].
This section is structured as follows: first, we describe the computational details used in common for the distorted hydrogen chain and iron crystals, including Hamiltonian construction, details of the quantum hardware employed and classical optimization schemes chosen.Then, we give details of the noise mitigation techniques applied for this work.Afterwards, we explain the system-dependent setups and simplifications that guide our choice of quantum algorithm for each model considered.

A. General Technical Details
Similarly to the molecular methodology, the PBCadapted second-quantized Hamiltonian requires a set of electronic integrals.We employ localized atomic Gaussian basis set adapted to translational symmetry [25].In contrast to the more popular use of plane wave basis sets in condensed matter simulations, it requires significantly smaller number of basis functions with the expensive integral evaluation costs as a drawback.Gaussian basis sets are favourable for quantum computations in the NISQ era, because a smaller number of basis functions leads to a smaller number of qubit and circuit depth requirements.We use the Jordan-Wigner scheme [32] to encode a chemistry Hamiltonian and an ansatz expressed with fermionic operators to the corresponding qubit operators to be used in Eq. ( 1) and Eq. ( 2), respectively [2,5].
The present calculations are dependent on a number of software packages.The periodic molecular integral evaluations associated with the Hartree-Fock (HF) calculations are performed by using the classical computational chemistry package PySCF version 1.7.2.[33].Gaussian density fitting is used for efficiently handling electron integrals.To facilitate convergence, primitive Gaussian functions with exponents less than 0.1 are discarded to avoid severe linear dependency.Divergence in exchange energy at the center of the Brillouin zone of PBCadapted HF or hybrid DFT is addressed by the Ewald correlation scheme.The same tool is used to perform reference (classical) post-HF calculations [34,35] such as k-point dependent Møller-Plesset perturbation theory (MP2) and coupled cluster singles and doubles (CCSD).PBC-adapted quantum computational methods (Tran-sQSE and PBC-adapted VQE) are implemented into our quantum computational chemistry software InQuanto [36], combined with the OpenFermion package version 0.11.0 [37] to facilitate fermionic operations.Quantum circuits are optimized and compiled for each quantum backend with the retargetable quantum compiler tket [38] via its Python interface, pytket version 0.6.0.Quantum simulations running on the classical computer are performed with Qulacs [39] for state-vector simulations and QasmSimulator of Qiskit [40] for shot-based sampling simulations, both via pytket interface.
To execute variational optimization, we employ two classical optimization algorithms: Rotosolve [41] and Stochastic Gradient Descent (SGD).[42,43].Rotosolve is a gradient-free optimization algorithm based on machine-learning techniques originally designed for hardware-efficient ansätze, whereas the SDG algorithm requires gradient evaluated also by the quantum computer [44].Variational optimization of TransQSE is performed by using Rotosolve, whereas VQE with UCCSD-PBC ansatz is performed by using both Rotosolve and SGD.
Quantum hardware experiments are performed by using ibmq casablanca via cloud service, which is one of the IBM Quantum Falcon Processors.The number of shots is chosen to be 24000 for each quantum circuit measurement process.
The approximate correlation energy E corr defined as the difference between post-HF and HF energies is the primary target in our quantum hardware experiments.For convenience, hereafter we use ∆E(θ) defined as where E • HF denotes the HF energy calculated with the classical computer and E total (θ) refers to the total energy calculated by using Eq. ( 1) with the chemistry Hamiltonian on quantum hardware or simulator, which may be influenced by the noise and/or stochastic error.If these errors are negligible, the value of ∆E(θ) with optimal parameters θ is equivalent to E corr , which must be a negative value.However, noise in the NISQ device can make it even a positive value, because noise-induced high energy excited states may contaminate the calculated ground state wavefunction.

B. Noise Mitigation Methods
By the nature of the chemistry-motivated ansätze [45] used in this study, the effect of hardware noise may be understood in terms of involvement of states that should not participate in the ground state wavefunction, like those not guaranteeing particle number conservation (ionized configurations) or involving a different spin sector ("spinflipped" states) than the one corresponding to the ground state.Noise mitigation techniques with these ansätze work to partially suppress these type of symmetry violations.
In the present paper we employ two types of noise mitigation techniques to post-process the resulting shot counts.The first technique is SPAM mitigation implemented in pytket, whereas the second technique is PMSV implemented in InQuanto.
SPAM mitigation [31] assumes that noise-induced errors are independent on the circuit to be executed, but only occur during the state preparation and measurement steps.Therefore the density matrix is preliminarily sampled to get the noise profile of the device.Then, application of the inverse of the density matrix to the quantum state would suppress the error caused by these noise channels.
PMSV is a novel technique we have developed to symmetry-verify quantum calculations in an efficient way.Molecular symmetries can be represented as a linear combination of Pauli strings.Symmetries such as mirror planes (Z 2 ) and electron-number conservation (U 1 ) can be represented by a single Pauli string that tracks the parity of the wavefunction [46].One can exploit these symmetries to apply qubit tapering techniques if the system is described by an Abelian point group [47,48].Alternatively, point group symmetry may be applied to mitigate noise.In PMSV, we perform symmetry verification by taking advantage of the commutation of the Pauli symmetries with terms in the Hamiltonian.We can partition these Hamiltonian terms and Pauli symmetries into groups of commuting Pauli strings [49,50].If each Pauli string in the partitioned set commutes with the symmetry operator, then each measurement circuit that measures that commuting set can be symmetry-verified without extra quantum resources, discarding the measurements that violate the point group symmetries of the system.Periodic systems are described by space groups, but by transforming into the reciprocal space, they can be described by point groups, so that we can apply PMSV for solid-state systems as in molecular systems [13].A more detailed explanation of this method can be found in Appendix A.
Compared to the other error mitigation methods, PMSV would be characterized as a method exploiting the symmetry with no extra quantum resources, which is scalable and independent of the type of noise.One can perform symmetry verification using mid-circuit measurements, ancilla readouts or the quantum subspace expansion algorithm [51], at the expense of extra quantum resources.Some studies using actual quantum computers [52,53] employ McWeeny purification for 2 qubit systems.The one-body reduced density matrix (1-RDM) can be easily purified, but the two-body reduced density matrix (2-RDM) cannot be extended to larger systems without further theoretical investigation.PMSV works with both 1-RDM and 2-RDM in a scalable manner.Other studies use a method relying on the extrapolation technique.For example, this study [54] assumes that noise in a basis of states is the same as that in the actual simulation, then applies an inversion matrix to clean the results.PMSV does not have any assumptions on the type of noise.

C. Setup for Hydrogen Molecular Lattices
Let us proceed to the first system to be investigated, a distorted hydrogen chain.It is usually taken as the first target of the PBC-adapted variational calculations as in Refs.[7,8,13,14,55].More extensive calculations using a quantum simulator running on the classical computer for systems including H 2 , He, and LiH lattices up to three-dimensional lattices are available in previous work [13].
We consider a one-dimensional hydrogen chain with alternating bond length with d and 1.5d, with d being equal to 0.75 Å (comparable to the molecular equilibrium H-H distance).With this distortion, we ensure the insulating character of the model.See Fig. 1 for the schematic representation of the geometry in the unit cell.
To generate the integrals in the crystalline orbitals, kpoint dependent restricted HF (RHF) with the STO-3G basis set is used.The spin multiplicity and the total charge per H 2 molecule are set to be 1 (spin-singlet) and 0 (charge-neutral), respectively.Here we consider a 2 k-point model.One could naively consider performing VQE hardware calculations with the UCCSD-PBC ansatz for the 2 kpoint model, requiring four qubits with six variational parameters to be optimized after applying the qubit tapering technique based on the symmetry of the system [47,48].However, this setup has been found not to be feasible on currently available hardware because of the large circuit depth of the corresponding ansatz, which includes many excitations with amplitudes of similar absolute value [13].
For this reason, instead of VQE with UCCSD-PBC ansatz, we demonstrate the TransQSE algorithm [13] with further approximations for this system.To perform the TransQSE experiment for the 2 k-point hydrogen chain, the Hamiltonian and the ansatz are transformed to a localized space representation spanned by Wannier orbitals [56].In this localized representation it has been found that a UCCSD-VQE wavefunction for the chosen distorted hydrogen chain geometry has two double excitations with significantly larger amplitudes than the others.This is in contrast with the momentum space solution, where no predominant excitations have been found.The dominance of these excitations can be also observed in the MP2 approximated amplitudes, and furthermore, the amplitudes of the dominating excitations are the same due to translational symmetry.By neglecting all other excitation operators apart from these two major ones, the cluster operator in the UCC ansatz in the real space ÛR = e T − T † is simplified as where â † R,p,σ and âR,p,σ are the creation and the annihilation operators for the localized orbital p with spin σ on the R-th H 2 molecule in the chain.Both excitations are inter-molecular and due to the periodic boundary the second term is obtained from the first term by applying a translation operator Λ, which shifts this term to the next unit cell on the right.By taking advantage of the translational operator, instead of the Trotterized UCC ansatz we write the TransQSE ansatz as where |Φ 0 and the ground-state energy is found by minimizing the energy function where as discussed in Ref. [13].Since the excitation only involves a limited range of 4 qubits in |Ψ W (θ) , the operators in Eq. ( 6) are reduced to only 4 qubits from 8 by contraction.Further resource reduction can be made by exploiting the symmetry [48].This system has a set of symmetry operators Ŝ as The first and second operators are associated with particle conservation for ↑ and ↓ spins, respectively, and the third one represents the mirror plane spatial symmetry [46].We use the first two operators for obtaining the 2qubit system by tapering off the third and fourth qubits, and apply the third operator to PMSV noise mitigation.After these simplifications are considered, one needs to measure the expectation value of five 2-qubit Pauli strings to determine E TQSE (θ) with the 2-qubit ansatz, U (θ)|Ψ 0 = e −iθ Ŷ0 X1 |0 0 0 1 .To further reduce the number of measurements, the commuting Pauli strings were measured together using a graph coloring algorithm implemented into pytket, totalling to two circuits to be measured.
To facilitate the variational optimization with hardware, the term 1/(1 + s 1 (θ)) in Eq. ( 6) is expanded into a Taylor series as As the overlap term S(θ) is expected to be small in magnitude due to the relatively low correlation, taking the first order of S(θ) should be a good approximation.It gives us a steeper function in the large |θ| region to accelerate the convergence.Hereafter we use Eq. ( 11) to the first order of S(θ) as the cost function of TransQSE.

D. Setup for Iron Crystals
To set a practical starting point of the quantum computational iron calculations, here we compute the energy difference between the body-centered cubic [BCC, Fig. 2(a)] and face-centered cubic [FCC, Fig. 2(b)] crystal structures, We focus on the ferromagnetic (FM) phase of iron in both BCC and FCC lattices, although it is known that the ground state of the FCC cell is antiferromagnetic (AFM) [18,20,23].This is a fair approximation because the energy difference between AFM-FCC and FM-FCC is fairly small [23].Thus, we consider only the FM-BCC structure with only one possible spin configuration, which is technically simpler as it allows us to consider a small simulation cell.Hereafter we refer to this BCC-FCC energy difference without mentioning their magnetic features.We shall leave the consideration of AFM-FCC structures with a variety of spin configurations for future studies.This is due to the increased resources required, as the model must be constructed with prohibitively large super-cells, and singledeterminant HF solutions are bad starting points for the calculations.
We exploit the simplified model system in the following two steps.i) Classical CCSD calculations are performed for the system with a relatively large number of k-points (which may be too many for VQE with currently available hardware) obtained to characterize the effect of electron correlation on the energy difference.ii) Simplified model systems with the smaller number of k-points and/or active orbitals are designed via a careful analysis of the larger system investigated in the step (i).In contrast to the hydrogen chain case, we identify a few leading excitations in the cluster operator in momentum space.This allows us to construct a simplified UCCSD-PBC ansatz for our VQE calculations.
The resulting simplified model system is subjected to the hardware experiment of VQE with UCCSD-PBC ansatz [13].The momentum-space second-quantized Hamiltonian ĤK is expressed as [13] ĤK = where h P Q and h P R QS are the transformed one-and twoelectron integrals of the periodic system, with P, Q, R, S denoting the composite indices of k-point in the Brillouin zone k, spatial orbitals p, and spin σ.
L3 b 3 with b 1 ,b 2 and b 3 denoting the reciprocal lattice vectors and k (1) , k (2) and k (3) are integers such that − Lα 2 < k (α) ≤ Lα 2 for α = 1, 2, 3 for the [L 1 L 2 L 3 ] k-point mesh.Index p labels orbitals in energy-ascending order of each k-point.The index is actually determined via a mapping function as P = q K (k, p, σ) [13].ĉ † P (ĉ P ) is the creation (annihilation) operator at the orbital P .The prime symbol in Eq. ( 12) on the sum indicates that it runs only with indices satisfying crystal momentum conservation [25].The cluster operator of the UCCSD-PBC ansatz ÛK = e T − T † is expressed as where t A I and t AB IJ are complex-valued coupled cluster amplitudes in general, with I, J and A, B being the composite indices for occupied and virtual orbitals, respectively.
DFT calculations are performed to support the CCSD calculations by checking that the properties are not significantly changed as the number of k-points increases.See Appendix B for the details.
We employ the basis set family of Los Alamos National Lab (LANL) effective core potentials (ECPs).[58] LANL-ECPs plus double-zeta basis set (LANL2DZ [59]) is used for all the iron calculations.Single-reference calculations are performed with the k-point dependent unrestricted HF (UHF) method implemented in PySCF.The number of unpaired electrons is set to be 2 per atom to represent the iron FM phases.Total charge is set to be neutral.
Technically, UHF calculations are not numerically stable.From our experience, this feature tends to be prominent in the larger lattice constant.If the initial guess of UHF is generated at each lattice constant, the potential energy curve may become discontinuous due to the different reference configuration.To obtain a smooth curve for the consistent analysis, we start our calculation with small lattice constant with a locally generated initial guess, and then use the converged crystal orbitals for calculating the next point in the lattice-constant-ascending order.

III. RESULTS AND DISCUSSION
A. TransQSE for distorted hydrogen molecular chain with 2 k-points

Variational optimization with quantum hardware
Here we demonstrate the variational optimization of TransQSE with the approximate cost function Eq. (11) to the first order by using Rotosolve with quantum hardware.Initial guess of the amplitude is set to be θ = 10 −5 (virtually HF state but with nonzero amplitude to avoid over simplification of the quantum circuit by tket).
The progress of the energy estimate is shown as lower triangles (connected by lines as visual guides) in Fig. although it should be virtually 0 kJ/mol by definition.Hereafter each uncertainty in the reported expectation values is calculated as the standard deviation evaluated from the samples for each θ taken from the noisy simulator emulating the hardware.PMSV can influence this uncertainty as some of the measurement outcomes may be discarded, but such an effect was found to be negligible in the present work, as the percentage of discarded shots is ∼ 1% for the H chain.The variational experiments converge in three iterations and find the optimal point ∆E(−0.0823)= −31 ± 4 kJ/mol, whose parameter θ is in good agreement with the state-vector result ∆E(−0.0928)= −46.03kJ/mol (the gradient is small d(∆E) dθ θ=−0.0928= 13 kJ/mol), but the energy is again shifted upwards due to the device noise.The error is larger than those of the iron model systems discussed in Sec.III B, probably because the energy is evaluated as a product of two expectation values as shown in Eq. ( 11) to enhance the relative error due to the noise.Thanks to SPAM correction, the resulting energy is driven to lower values than the classically evaluated HF energy, that is, ∆E < 0. The SPAM correction improves the total energy in all the points in Fig. 3 by about 50 kJ/mol on average.PMSV is found to display no significant effect on energy in this particular case.This is due to the relatively small value of θ, and SPAM already mitigating the effect of the symmetry-violating measurement results to be excluded by PMSV.See Appendix A if one is interested in the details.
As shown above, even in this simple case, it was found to be difficult to obtain a quantitative agreement between the energies obtained from state-vector simulations and the quantum device results used in these calculations.However, the optimal θ is reproduced accurately within 12% of relative error.Improvements in the agreement between theoretical and experimental ∆E(θ) mainly depend on quantum hardware development, although advances in noise mitigation techniques will also have an impact on these results.In contrast to the PBC-adapted VQE, TransQSE can still be compact even with more kpoints by truncating spatial long-range correlations.The state-vector simulation of TransQSE demonstrates that good approximations of the ground state energy value are obtained (∆E = −46.03kJ/mol) with a smaller number of parameters than the analogous PBC-adapted VQE calculation (∆E = −57.73kJ/mol [13]).However, exactly predicting in which cases TransQSE is a more feasible choice is generally not trivial.In this sense, cheap preprocessing techniques to screen important angles in the ansatz may help to solve this question.

Classical CCSD calculation with a [4 4 4] k-point mesh
First we investigate the effect of electronic correlation on the BCC-FCC energy difference by performing classical UHF/CCSD calculations.We will consider a model for these iron phases simplified in such a way that the most important features in these relatively complex systems are captured with as much fidelity as possible.
Within the available computational resources, we consider a UHF/CCSD setup with a [4 4 4] Monkhorst-Pack mesh, with the active space consisting of the crystal orbitals with orbital and spin indices (p, σ) = {(8, ↑), (9, ↑ ), (6, ↓), (7, ↓)} for each k value.These crystal orbitals correspond to linear combinations of iron d orbitals in the region around the Fermi level.The total number of active electrons and crystal orbitals are equal to 128 and 256, respectively.The reference electronic configuration is not uniform between k-points, but each k sector contains 0 to 4 electrons in it.We have indirectly confirmed that this [4 4 4] k-point mesh is sufficient by using DFT with the same basis set, because the potential energy curves were not significantly changed as the number of k-points increased (see Appendix B).The number of k-points is still significantly smaller than the converged values in literature [19,23,60,61], but it is large enough for the k-point sampling error to be sufficiently small for the purposes of this study.
One of the remarkable features in the comparison of UHF vs. CCSD results is that UHF significantly overestimates the BCC-FCC energy difference, but CCSD corrects this energy difference to a value comparable to accurate plane-wave basis DFT results [23], which are used as reference values.This trend is illustrated in our calculations of the potential energy curve as a function of the lattice constant with UHF for each of FCC and BCC, as shown in Fig. 4. The UHF calculations accurately reproduce the reference lattice constants, but the estimated BCC-FCC difference (30 kJ/mol) is significantly overestimated, to almost twice the reference value as shown in Fig. 4.Such a disagreement is corrected by taking the electronic correlation into account with CCSD.The total energy with CCSD/UHF at the optimized lattice constant with UHF for each of BCC and FCC is shown in Fig. 4. We find that the energy difference is corrected to 12 kJ/mol, which is indeed comparable to the reference value indicated by the horizontal line in Fig. 4.
These results suggest that the FCC phase is more influenced by the electronic correlation than the BCC phase.We will not provide a detailed interpretation, but many competing phases may contribute to increasing the electronic correlation in the FCC case, in contrast to the isolated BCC phase [19].Hereafter we focus on this feature to construct simplified model systems to set a realistic starting point for experiments on currently available quantum computers.In other words, our simplifications should take into account that the correlation energy in our simplified FCC model will be larger than that of the BCC model.

Simplified system with [2 1 1] k-point mesh
As a first step towards the quantum hardware experiments, we start from a simplified model system to adjust the resources needed for PBC-adapted VQE calculations to the constraints imposed by the device we are going to target.Let us employ the simplest possible [2 1 1] k-point mesh to test if this simplified system can reproduce the FCC-BCC features mentioned in Sec.III B 1.
The potential energy curves calculated with UHF on this model [Fig.5] indicate that UHF overestimates the BCC-FCC energy difference, which is the same trend observed in the [4 4 4] k-point mesh cases.Although these simplified models are quantitatively insufficient, the UHF calculation returns a good estimate of the equilibrium lattice constant.This implies that this drastically simplified model is able to qualitatively reproduce the properties of the iron phases in which we are interested.
The potential energy curves calculated with UHF/CCSD, including all the excitation operators in the cluster operator T , are shown in Fig. 5.Although the BCC-FCC energy difference is still overestimated, it is indeed corrected in the right direction.As shown in Fig. 5, similar results are obtained by using the same active space used for the [4 4 4] k-point mesh case, that is, (p, σ) = {(8, ↑), (9, ↑), (6, ↓), (7, ↓)} for each k.The reference configuration is given as {1, 1, 1, 1} for k = 0 and {0, 0, 0, 0} for k = 1, where k ≡ k = k 2 b 1 .This suggests that the correlation energy is dominated by the frontier orbitals as in the molecular cases [62].Therefore, TABLE I: Approximate correlation energy Ecorr of iron crystal calculations with 2 k-points obtained by using either CASCI, CCSD or VQE with the UCCSD-PBC ansatz.The active space used is identified: {(k, 8, ↑), (k, 9, ↑), (k, 6, ↓), (k, 7, ↓)} for k = 0, 1 (see Fig. 6).this simple model is validated and we will employ the [2 1 1] k-point mesh with the active space for quantum hardware experiments.Hereafter the optimized lattice constants (2.84 Å for BCC and 3.64 Å for FCC) obtained by UHF/CCSD with the active space are used.The resulting approximate correlation energy E corr for each of BCC and FCC is shown in Tab.I, which is used for comparing with the quantum computational methods.

Method excitation operator(s) Ecorr
It is required to further simplify the model system, because the quantum circuits corresponding to the PBCadapted VQE with the UCCSD ansatz are still too deep for the currently available quantum hardware.To demonstrate the quantum hardware calculations, let us perform further simplification of the system based on the same idea applied for the TransQSE calculations, that is, focusing on the excitation operators with prominent amplitudes.We performed state-vector simulations of VQE for iron with the active space mentioned above.The resulting energies of both BCC and FCC are in good agreement with the CCSD counterpart (see Tab. I).The complete active space configuration interaction (CASCI) via exact diagonalization is also performed for comparison.The nonzero cluster amplitudes of VQE are listed in Tab.II and Tab.III for BCC and FCC, respectively.Note that all the amplitudes in the ansatz are of real number as a consequence of the reduction to 2 k-points.In both cases, we can find one prominent cluster amplitude corresponding to a inter-k-point double excitation.This inter-k-point double-excitation has more than ten times larger amplitude in magnitude compared to the others.Therefore we may define an even more compact active space consisting of these four crystal orbitals with an ansatz including one double excitation.
As shown in Tab.I, more than 90% of the reference correlation energy is reproduced with this one double excitation operator in both BCC and FCC cases.More excitation operators in the ansatz can in principle improve the correlation energy, but such a small difference may not be resolved with the deepened quantum circuit on the device currently available.Therefore the cluster operator is approximately expressed as TABLE II: Non-zero cluster amplitudes t AB IJ of the PBCadapted VQE calculation with the UCCSD ansatz for the 2 k-point model of iron BCC structure with the active space (see text).As a consequence, we have elaborated a simple computational setup consisting of four active crystal orbitals with one double-excitation operator of the UCCSD-PBC ansatz with two electrons.See Fig. 6 for the schematic representation of the orbital energy levels and the active space.The same symmetry operators as in the model system for TransQSE (Eq.( 10)) are used for reducing the sytem to an equivalent 2-qubit system by applying the qubit tapering technique, and for PMSV noise mitigation.To determine the expectation value Ψ(θ)| ĤK |Ψ(θ) , one needs to measure the expectation value of four 2-qubit Pauli strings with the ansatz U (θ)|Ψ 0 = e −iθ Ŷ0 X1 |0 0 0 1 .Applying the same partitioning by using pytket, two circuits are actually measured in both BCC and FCC models.

Variational quantum hardware experiments
Let us proceed to demonstrate the variational optimization of the energy by using real quantum hardware provided through cloud by IBM Quantum.We need at least ∼10 kJ/mol of accuracy to reliably discuss the BCC-FCC energy difference, but reaching this threshold has been found to be unfeasible without noise mitigation in our preliminary calculations, even for this simplified 2-qubit model system.To realize this accuracy on the NISQ device, we apply SPAM and PMSV noise mitigation.By using these techniques for ibmq casablanca, a sufficiently accurate expectation value is estimated within ∼ 5 kJ/mol of error compared to the state-vector simulation results.
First we show the results of variational optimization with Rotosolve.The progress of correlation energy along the variational optimization by using Rotosolve is shown in Fig. 7(a) and Fig. 7(b) for BCC and FCC, respectively.In both cases, optimization rapidly converges to the point that it is in very good agreement with the energy obtained by the state-vector simulation.The resulting parameter and the energy are represented as ∆E(−0.5561)= −215 ± 3 kJ/mol for BCC and ∆E(−0.6522)= −265 ± 2 kJ/mol for FCC.The small number of optimization steps owes to the properties of Rotosolve, which is designed for systems described by quantum circuits with parameterized rotation gates.reliably tell the difference in correlation energy (not the total energy) between FM-BCC and FM-FCC (51 kJ/mol), in which FM-FCC has more correlation energy to correct the UHF value indeed to the right direction.
Similarly, the SDG optimizer can also accurately locate the ground state obtained in the reference statevector simulation.The results of the progress of ∆E(θ) along the variational optimization are shown in Fig. 8(a) for BCC and Fig. 8(b) for FCC.The resulting energies are also in good agreement with the exact results from state-vector simulations.In contrast to Rotosolve, the SDG method uses energy derivatives with respect to the circuit parameter, which is also evaluated by using the quantum computer.The resulting parameter and energy are represented as ∆E(−0.5196)= −214 ± 3 kJ/mol for BCC and ∆E(−0.6060)= −264 ± 2 kJ/mol for FCC.The resulting energy of SDG agrees with that of Rotosolve for each of BCC and FCC.The parameters look different, but it is likely due to the small gradient around the optimal point, which is implied by the slow improvement of the energy in Fig. 8.The SDG optimizer takes more steps in this 1-parameter case, but it may become more efficient if we take more excitation operators (circuit parameters) into account in the ansatz.As in the case of the hydrogen chain, the percentage of discarded shots due to the PMSV correction is low, of the order of ∼ 5% for this system.
In this particular case, additional excitation operators to the ansatz would not have a significant effect on the energy (∼5% of correlation improvement compared to the first one).Larger hardware experiments with more qubits and more excitation operators (i.e., wider and deeper quantum circuits) with further noise mitigation techniques including PMSV are the subject of future work, and will be reported elsewhere.

IV. CONCLUSIONS AND OUTLOOK
We have performed quantum hardware calculations for the solid-state model systems with periodic boundary conditions (PBC).We employed a distorted hydrogen chain and iron crystals as the targets, which are both systematically simplified to set a starting point for simulations on current quantum hardware with 2 qubit 1 parameter ansätze.We have applied the TransQSE method for the hydrogen chain and the PBC-adapted VQE method for the iron crystals.
To make the most of NISQ devices, we have applied noise mitigation techniques to improve the agreement between the experimental energy estimate and reference state-vector simulation values.We have applied the state preparation and measurement (SPAM) correction for all the variational optimizations on hardware.In all the cases the energy estimates were improved.In addition to SPAM correction, we applied the novel PMSV noise mitigation technique, in which shot counts are post-selected so that the intrinsic symmetries of the system (Z 2 and U 1 ) are not violated.SPAM and PMSV are simultaneously used to realize ∼ 5 kJ/mol of agreement in the calculations of simplified iron crystal models.Although these results are for very simple model systems, we be- lieve that these results set an important starting point for systematic improvement of quantum chemical calculations on quantum computers by rolling back the simplification procedure presented in this paper.Once quantum computers improve, larger basis sets and k-point grids can be considered, and much more accurate estimates of the total energy for these systems can be obtained.
In addition to the hardware improvement itself, several techniques would be useful to extend the range of periodic systems that we may simulate with quantum computers in the near term.One of the promising approaches to handle a large system is using quantum embedding techniques like density matrix embedding theory (DMET) [63] or Dynamical Mean-Field Theory [64,65].For example, the complete unit cell energy in DMET is calculated as the sum of fragment energies.The interaction between each fragment and the rest of the unit cell (referred to as bath) is calculated self-consistently.The flexibility of the choice of the size of the fragment allows us to reduce the quantum resource requirements to address bigger unit cells step by step.Similar gains can be obtained by using other embedding techniques.
We highlight that consideration of the translational nature of periodic systems, both in direct and reciprocal space, has allowed us to apply simplifications that mitigate the additional resources needed to treat our extended systems on a quantum computer.This is a significant advantage over the consideration of molecular systems with similar sizes in terms of basis sets and numbers of qubits required to perform a simulation.More gen-erally, exploitation of quantum resource savings derived from translational symmetry should always be taken into account when one deals with PBC systems.Now that all the basic tools to perform PBC-adapted VQE calculations have been tested on hardware, we expect that larger numbers of basis functions and k-points with wider active spaces will be possible as hardware improves.Among others, bigger hardware calculations with more complex ansätze and a variety of noise mitigation techniques are underway in order to go one step ahead, which will be reported elsewhere.
Symmetry verification in the NISQ era aims to postselect on quantum measurement results based on a measurable property of a quantum circuit during runtime of the computation.The verification procedure reads out a conserved property during the quantum calculation with respect to the projection P expressed as where Î is the identity operation, and Ŝ is a symmetry operation with the parity represented by (−1) x , with x being an integer.Techniques in literature include symmetry verification via mid-circuit measurements or ancillacontrolled operations plus readouts [51].The inherent problem with these techniques is the introduction of additional two-qubit gate operations.Another difficulty is identifying what symmetries a system possesses and how to represent them as operations on many-qubit systems, which may be a non-trivial process.For applications in quantum chemistry, point-group symmetries and spacegroup symmetries can be used to define the symmetries for systems of interest.These can be represented as operations on many-qubit systems by following prescribed routines available in literature [48].
We introduce a technique that can symmetry-verify a quantum circuit without imposing the extra two-qubit gate penalty.In the NISQ era, this is a necessity for any simulation that aims to preserve particular properties.
An example would be the use of quantum simulations to find ground states for molecular and material science problems.Our symmetry verification routine makes use of a measurement reduction technique implemented in tket.This technique finds partitions of commuting Pauli strings and maps them to measurement circuits [50].This method is referred to as Partition-Measurement Symmetry Verification (PMSV).
In PMSV, one partition of commuting Pauli-strings corresponds to one measurement circuit.To make a measurement circuit "symmetry-verifiable", we check if the partition of commuting Pauli strings commutes with every Pauli symmetry that a chemical system exhibits.If this check is true, these Pauli-symmetries can be measured by the measurement circuit as well as the Pauli strings of the Hamiltonian.When we process the quantum measurement result, we measure the parity, x, of our Pauli-symmetry first.We use this check as a constraint to keep or discard a particular measurement result, before post-processing the expectation value of the other Pauli strings.We observe that our symmetry verification scheme can be performed for any quantum simulation via the application of a small classical processing step, but without extra quantum requirements.
The algorithm to prepare a set of quantum circuit and result mapping is given in Algorithm 1, whereas that to post-select the measurement result is shown in Algorithm 2.

Algorithm 1 Prepare Symmetry Verifiable
Measurement Circuits We demonstrate how PMSV works by comparing the probabilities of the measurement outcome from the noisy simulation with different error mitigation techniques.Single point energy calculation with the preliminary optimized parameter (by state-vector simulation) is considered for convenience.Noise model is taken from the IBM

FIG. 1 :
FIG. 1: Schematic representation of the unit cell (indicated by dashed lines) of the distorted hydrogen chain with alternating bond lengths, d and 1.5d.The bond parameter is set to be d = 0.75 Å in the present paper.

3 .FIG. 3 :
FIG. 3:Progress of energy estimate ∆E(θ) and along the variational optimization on hardware (ibmq casablanca) by using the rotosolve algorithm (solid lines with markers).(E total , E • HF ) = (−2.7226×10 3 kJ/mol, −2.6766×10 3 kJ/mol) according to the state-vector simulations.The dashed line denotes the correlation energy Ecorr for the given model Hamiltonian.The optimization is performed for the expectation values calculated with SPAM-and PMSV-corrected shot counts (lower triangles).Those calculated with raw (circles) and SPAM-corrected (upper triangles) shot counts are also shown for comparison.

FIG. 4 :
FIG.4: Potential energy curves of FM-BCC and FM-FCC iron as functions of lattice constants for[4 4 4]  k-point mesh at UHF level.Also represented are single-point calculations at CCSD level at optimized lattice constant values for BCC and FCC.The energy baseline is set so that the lowest energy is equal to zero for each UHF and UHF/CCSD.The optimized lattice constants from the accurate plane-wave basis DFT calculations[23] are indicated by vertical lines (BCC: 2.83 Å, FCC: 3.64 Å), while the BCC-FCC energy difference is represented by the horizontal red line (14.67 kJ/mol).

FIG. 5 :
FIG. 5: Potential energy curves as functions of lattice constants for [2 1 1] k-point mesh with UHF, CCSD/UHF with active space [indicated by CCSD(a), See Fig. 6], and CCSD/UHF with all the valence orbitals [indicated by CCSD(b)].UHF overestimates the BCC-FCC energy difference (168 kJ/mol), and CCSD corrects it to some extent [118 kJ/mol for CCSD(a) and 100 kJ/mol for CCSD(b)], which is qualitatively the same feature found in the accurate case with [4 4 4] k-points.The meanings of the vertical and horizontal lines in each panel are the same as those in Fig. 4.

7 ,FIG. 6 :
FIG. 6: Schematic representation of the orbital energy levels for the simplified 2 k-point iron model.Each horizontal bar represents the energy level with the number corresponding to the orbital index p with spin σ = {↑, ↓}.Occupied orbitals are indicated by the points "•".The active space used in the present hardware calculations is represented by the red rectangle in color.

FIG. 7 :
FIG. 7: Evolution of the ∆E with respect to the number of steps used by the Rotosolve optimization process.(E total , E • HF ) = (−3.228573× 10 5 kJ/mol, −3.226368 × 10 5 kJ/mol) for BCC and (E total , E • HF ) = (−3.227450× 10 5 kJ/mol, −3.224735 × 10 5 kJ/mol) for FCC according to the state-vector simulations.Black points show hardware results obtained with the ibmq casablanca device.The line is included as a visual guide to show the progress of the correlation energy estimate along the variational optimization.The dashed line is the correlation energy Ecorr for the given model Hamiltonian.FCC result has 51 kJ/mol more correlation energy in magnitude to correct the BCC-FCC energy difference from 168 kJ/mol to 117 kJ/mol, which is qualitatively the same feature found in the accurate case with [4 4 4] k-point mesh.SPAM and PMSV are applied at the same time.Raw and SPAM-corrected ∆E are given for each optimization step to show the effect of noise mitigation.

FIG. 8 :
FIG. 8: Evolution of ∆E with respect to the number of steps used by the SDG optimization process.(E total , E • HF ) = (−3.228573× 10 5 kJ/mol, −3.226368 × 10 5 kJ/mol) for BCC and (E total , E • HF ) = (−3.227450× 10 5 kJ/mol, −3.224735 × 10 5 kJ/mol) according to the state-vector simulations.Black points show hardware results obtained with the ibmq casablanca device.The black line is included as a visual guide to show the progress of the correlation energy estimate along the variational optimization.The dashed line is the correlation energy Ecorr for the given model Hamiltonian.Error mitigation based on SPAM and PMSV are applied at the same time.Raw and SPAM-corrected ∆E are given for each optimization step to show the effect of noise mitigation.

TABLE III :
Non-zero cluster amplitudes t AB IJ of the PBCadapted VQE calculation with the UCCSD ansatz for the 2 k-point model of iron FCC structure with the active space (see text).