SQD para estimativa de energia de um hamiltoniano químico

Nesta aula, aplicaremos a SQD para estimar a energia do estado fundamental de uma molécula.

Em particular, discutiremos os seguintes tópicos usando a abordagem do padrão Qiskit $4$ -step:

1. Etapa 1: Mapear o problema para circuitos e operadores quânticos
   • Configure o Hamiltoniano molecular para $N_2$.
   • Explique o cluster unitário local Jastrow (LUCJ) inspirado na química e amigável ao hardware [1]
2. Etapa 2: otimizar para o hardware de destino
   • Otimizar a contagem de portas e o layout do ansatz para execução em hardware
3. Etapa 3: Executar no hardware de destino
   • Execute o circuito otimizado em uma QPU real para gerar amostras do subespaço.
4. Etapa 4: Resultados do pós-processamento
   • Introduzir o loop de recuperação de configuração autoconsistente [2]
     • Pós-processar o conjunto completo de amostras de bitstring, usando o conhecimento prévio do número de partículas e a ocupação orbital média calculada na iteração mais recente.
     • Criar probabilisticamente lotes de subamostras a partir de cadeias de bits recuperadas.
     • Projete e diagonalize o Hamiltoniano molecular em cada subespaço amostrado.
     • Salve a energia mínima do estado fundamental encontrada em todos os lotes e atualize a ocupação orbital média.

Usaremos vários pacotes de software durante a aula.

• PySCF para definir a molécula e configurar o Hamiltoniano.
• ffsim para construir o ansatz LUCJ.
• Qiskit para transpilar o ansatz para execução em hardware.
• Qiskit IBM Runtime para executar o circuito em uma QPU e coletar amostras.
• Qiskit addon SQD recuperação de configuração e estimativa de energia do estado fundamental usando projeção de subespaço e diagonalização de matriz.

1. Mapear o problema para circuitos e operadores quânticos

Hamiltoniano molecular

Um Hamiltoniano molecular assume a forma genérica:

$\hat{a}^\dagger_{p\sigma}$ / $\hat{a}_{p\sigma}$ são os operadores de criação/aniquilação fermiônica associados ao $p$ -ésimo elemento do conjunto de base e ao spin $\sigma$. $h_{pr}$ e $(pr|qs)$ são as integrais eletrônicas de um e dois corpos. Usando o site pySCF,, definiremos a molécula e calcularemos as integrais de um e dois corpos do Hamiltoniano para o conjunto de bases 6-31g.

import warnings
import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],  # Two N atoms 1 angstrom apart
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)  # hcore: one-body integrals
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)  # eri: two-body integrals

# Compute exact energy for comparison
exact_energy = cas.run().e_tot

Output:

converged SCF energy = -108.835236570774
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

Nesta lição, usaremos a transformação de Jordan-Wigner (JW) para mapear uma função de onda fermiônica para uma função de onda qubit, de modo que ela possa ser preparada usando um circuito quântico. A transformação JW mapeia o espaço de Fock dos férmions em M orbitais espaciais no espaço de Hilbert de 2M qubits, ou seja, um orbital espacial é dividido em dois orbitais de spin, um associado a um elétron com spin up ( $\alpha$ ) e outro com spin down ( $\beta$ ). Um orbital de spin pode estar ocupado ou desocupado. Normalmente, quando nos referimos ao número de orbitais, estamos usando o número de orbitais espaciais. O número de orbitais de spin será o dobro. Nos circuitos quânticos, representaremos cada orbital de spin com um qubit. Assim, um conjunto de qubits representará spin-up ou $\alpha$ -orbitals, e outro conjunto representará spin-down ou $\beta$ -orbitals. Por exemplo, a molécula $N_2$ para o conjunto de base 6-31g tem orbitais espaciais $16$ (ou seja, $16$ $\alpha$ + $16$ $\beta$ = $32$ orbitais de spin). Portanto, precisaremos de um circuito quântico de $32$ -qubit (podemos precisar de qubits de ancilla extras, conforme discutido posteriormente). Os qubits são medidos em bases computacionais para gerar cadeias de bits, que representam configurações eletrônicas ou determinantes (Slater). Ao longo desta lição, usaremos os termos bitstrings, configurações e determinantes de forma intercambiável. As cadeias de bits nos informam a ocupação de elétrons em orbitais de spin: um $1$ em uma posição de bit significa que o orbital de spin correspondente está ocupado, enquanto um $0$ significa que o orbital de spin está vazio. Como os problemas de estrutura eletrônica preservam as partículas, apenas um número fixo de orbitais de spin deve ser ocupado. A molécula $N_2$ tem elétrons $5$ spin-up ( $\alpha$ ) e $5$ spin-down ( $\beta$ ). Portanto, qualquer bitstring que represente os orbitais $\alpha$ e $\beta$ deve ter cinco $1\text{s}$ cada para a molécula $N_2$.

1.1 Circuito quântico para geração de amostras: a abordagem LUCJ

Nesta lição, usaremos o ansatz de Jastrow (LUCJ) \ [1] de cluster acoplado unitário local para a preparação do estado quântico e a amostragem subsequente. Primeiro, explicaremos os diferentes blocos de construção do ansatz UCJ completo e as aproximações feitas na versão local dele. Em seguida, usando o pacote ffsim, construiremos o ansatz LUCJ e o otimizaremos usando o transpilador Qiskit para execução em hardware.

O ansatz UCJ tem a seguinte forma (para um produto de $L$ camadas ou repetições do operador UCJ)

onde $\vert \Phi_{0} \rangle$ é um estado de referência, geralmente considerado o estado Hartree-Fock (HF). Como o estado Hartree-Fock é definido como tendo os orbitais de menor número ocupados, a preparação do estado HF envolverá a aplicação de portas X para definir os qubits correspondentes aos orbitais ocupados como um. Por exemplo, o bloco de preparação do estado HF para 4 orbitais espaciais e 2 de spin para cima e 2 de spin para baixo pode ter a seguinte aparência:

Uma única repetição do operador UCJ ${(e^{K^{(\mu)}} \times {e^{iJ^{(\mu)}}} \times {e^{-K^{(\mu)}}})}$ consiste em uma evolução de Coulomb diagonal ( $e^{iJ^{(\mu)}}$ ) intercalada por rotações orbitais ( $e^{K^{(\mu)}}$ e $e^{-K^{(\mu)}}$ ).

Os blocos de rotação orbital funcionam em uma única espécie de spin ( $\alpha$ (up-spin)/ $\beta$ (down-spin)). Para cada espécie de elétron, a rotação orbital consiste em uma camada de portas de um único qubit $R_{z}$ seguida por uma sequência de portas de rotação de dois qubits de Given (portas $XX + YY$ ).

As portas de 2 qubits atuam em orbitais de spin adjacentes (qubits vizinhos mais próximos) e, portanto, são implementáveis em IBM® QPUs sem a necessidade de portas SWAP.

O $e^{iJ^{(\mu)}}$, também conhecido como operador de Coulomb diagonal, consiste em três blocos. Dois deles funcionam nos mesmos setores de spin ( $e^{iJ_{\alpha \alpha}^{(\mu)}}$ e $e^{iJ_{\beta \beta}^{(\mu)}}$ ) e um funciona entre dois setores de spin ( $e^{iJ_{\alpha \beta}^{(\mu)}}$ ).

Todos os blocos em $e^{iJ^{(\mu)}}$ consistem em portas de número-número $U_{nn}(\phi)$ [1]. Uma porta $U_{nn}(\phi)$ pode ser dividida em uma porta $R_{ZZ}(\frac{\phi}{2})$ seguida por duas portas $Rz(-\frac{\phi}{2})$ de um único qubit que atuam em dois qubits separados.

Os componentes de mesmo spin ( $J_{\alpha \alpha}$ e $J_{\beta \beta}$ ) têm portas $U_{nn}$ entre todos os pares possíveis de qubits. No entanto, como as QPUs supercondutoras têm conectividade restritiva, os qubits devem ser trocados para realizar portas entre qubits não adjacentes.

Por exemplo, considere o seguinte bloco $e^{iJ_{\alpha \alpha}^{(\mu)}}$ (ou $e^{iJ_{\beta \beta}^{(\mu)}}$ ) para $N = 4$ orbitais espaciais. Para uma conectividade de qubit linear, as três últimas portas não são diretamente implementáveis, pois funcionam entre qubits não adjacentes (por exemplo, Q0 e Q2 não estão diretamente conectados). Portanto, precisamos de portas SWAP para torná-las adjacentes (a figura a seguir mostra um exemplo com $3$ portas SWAP).

Em seguida, o $J_{\alpha \beta}$ implementa portas entre os mesmos orbitais indexados de diferentes setores de spin (por exemplo, entre $0\alpha$ e $0\beta$ ). Da mesma forma, se os qubits não forem fisicamente adjacentes em uma QPU, essas portas também exigirão SWAPs.

Com base na discussão acima, o ansatz UCJ enfrenta alguns obstáculos para a execução de HW, pois precisa de portas SWAP devido a interações de qubit não adjacentes. A variante local do ansatz UCJ, LUCJ, aborda esse desafio removendo alguns $U_{nn}$ do operador de Coulomb diagonal.

Nos mesmos blocos de espécies de elétrons, $J_{\alpha \alpha}$ e $J_{\beta \beta}$ ), mantemos apenas as portas $U_{nn}$ compatíveis com a conectividade do vizinho mais próximo e removemos as portas entre qubits não adjacentes na versão LUCJ. A figura a seguir mostra o bloco LUCJ após a remoção de portas não adjacentes.

Em seguida, a versão LUCJ do bloco $J_{\alpha \beta}$ que funciona entre diferentes espécies de elétrons pode assumir formas diferentes com base na topologia do dispositivo.

Aqui, também, a versão LUCJ elimina as portas não compatíveis. A figura abaixo mostra variantes do bloco $J_{\alpha \beta}$ para diferentes topologias de qubit, incluindo grade, hexagonal, heavy-hex e linear.

• Quadrado : podemos ter portas $U_{nn}$ entre todos os orbitais $\alpha$ e $\beta$ sem nenhum SWAP e, portanto, não precisamos remover nenhuma porta $U_{nn}$.
• Heavy-hex : As interações $\alpha$ - $\beta$ são mantidas entre cada $4$ -ésimo orbital de spin indexado (como o 0º, 4º e 8º) e são mediadas por ancilla, ou seja, precisamos de qubits de ancilla entre as cadeias lineares que representam os orbitais $\alpha$ e $\beta$. Esse arranjo precisa de um número limitado de SWAPs.
• Hexagonal : Todos os outros orbitais, como o 0º, 2º e 4º orbitais indexados, tornam-se vizinhos mais próximos quando $\alpha$ e $\beta$ são dispostos em duas cadeias lineares adjacentes.
• Linear : Apenas um orbital $\alpha$ e um $\beta$ estão conectados, o que significa que o bloco $J_{\alpha \beta}$ terá apenas uma porta.

Diagramas de

Embora a remoção de portas da ansatz UCJ para construir a versão LUCJ a torne mais compatível com HW, a ansatz perde um pouco de expressividade. Portanto, mais repetições ( $L$ ) do operador UCJ modificado podem ser necessárias ao usar o ansatz LUCJ.

1.2 Inicialização do LUCJ ansatz

O LUCJ é um ansatz parametrizado, e precisamos inicializar os parâmetros antes da execução do hardware. Uma maneira de inicializar o ansatz é usar as amplitudes t1 e t2 do método CCSD (cluster singles and doubles) clássico acoplado, em que as amplitudes t1 são o coeficiente de operadores de excitação simples e as amplitudes t2 são para operadores de excitação dupla.

Observe que, embora a inicialização do ansatz LUCJ com t1 e t2 amplitudes gere resultados decentes, os parâmetros do ansatz podem precisar de mais otimização.

# Get CCSD t2 amplitudes for initializing the ansatz
ccsd = pyscf.cc.CCSD(
    scf, frozen=[i for i in range(mol.nao_nr()) if i not in active_space]
)
ccsd.run()

t1 = ccsd.t1
t2 = ccsd.t2

Output:

E(CCSD) = -109.0398256929733  E_corr = -0.20458912219883

1.3 Construindo a abordagem LUCJ usando ffsim

Usaremos o pacote ffsim para criar e inicializar o ansatz com t1 e t2 amplitudes calculadas acima. Como nossa molécula tem um estado Hartree-Fock de casca fechada, usaremos a variante de equilíbrio de spin do ansatz UCJ, UCJOpSpinBalanced.

Como o hardware do IBM tem uma topologia heavy-hex, adotaremos o padrão em zigue-zague usado em [1] e explicado acima para interações de qubit. Nesse padrão, os orbitais (qubits) com o mesmo spin são conectados com uma topologia de linha (círculos vermelhos e azuis). Devido à topologia heavy-hex, os orbitais de diferentes spins têm conexões entre cada quarto orbital, ou seja, o 0º, o 4º, o 8º e assim por diante (círculos roxos).

import ffsim
from qiskit import QuantumCircuit, QuantumRegister

n_reps = 2
alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]
alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]

ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(
    t2=t2,
    t1=t1,
    n_reps=n_reps,
    interaction_pairs=(alpha_alpha_indices, alpha_beta_indices),
)

nelec = (num_elec_a, num_elec_b)

# create an empty quantum circuit
qubits = QuantumRegister(2 * num_orbitals, name="q")
circuit = QuantumCircuit(qubits)

# prepare Hartree-Fock state as the reference state and append it to the quantum circuit
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)

# apply the UCJ operator to the reference state
circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)
circuit.measure_all()
# circuit.decompose().draw("mpl", scale=0.5, fold=-1)

O ansatz LUCJ com camadas repetidas pode ser otimizado com a fusão de alguns blocos adjacentes. Considere um caso para n_reps=2. Os dois blocos de rotação orbital no meio podem ser fundidos em um único bloco de rotação orbital. O pacote ffsim tem um gerenciador de passagem chamado ffsim.qiskit.PRE_INIT para otimizar o circuito mesclando esses blocos adjacentes.

2. Otimize para o hardware de destino

Primeiro, buscamos um backend de nossa escolha. Otimizaremos nosso circuito para o backend e, em seguida, executaremos o circuito otimizado no mesmo backend para gerar amostras para o subespaço.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService()
# Use the least-busy backend or specify a quantum computer using the syntax commented out below.
backend = service.least_busy(operational=True, simulator=False)
# backend = service.backend("ibm_brisbane")

Em seguida, recomendamos as seguintes etapas para otimizar o ansatz e torná-lo compatível com o hardware.

• Selecione os qubits físicos (initial_layout) do hardware de destino que adere ao padrão em ziguezague (duas cadeias lineares com um qubit de ancila entre elas) descrito acima. A disposição dos qubits nesse padrão resulta em um circuito compatível com hardware eficiente com menos portas.
• Gere um gerenciador de passes em etapas usando a função generate_preset_pass_manager do Qiskit com sua escolha de backend e initial_layout.
• Defina o estágio pre_init de seu gerenciador de passes em etapas para ffsim.qiskit.PRE_INIT. ffsim.qiskit.PRE_INIT inclui passagens do transpilador Qiskit que decompõem as portas em rotações orbitais e, em seguida, mesclam as rotações orbitais, resultando em menos portas no circuito final.
• Execute o gerenciador de passes em seu circuito.

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

spin_a_layout = [0, 14, 18, 19, 20, 33, 39, 40, 41, 53, 60, 61, 62, 72, 81, 82]
spin_b_layout = [2, 3, 4, 15, 22, 23, 24, 34, 43, 44, 45, 54, 64, 65, 66, 73]

initial_layout = spin_a_layout + spin_b_layout

pass_manager = generate_preset_pass_manager(
    optimization_level=3, backend=backend, initial_layout=initial_layout
)

# without PRE_INIT passes
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/o pre-init passes): {isa_circuit.count_ops()}")

# with PRE_INIT passes
# We will use the circuit generated by this pass manager for hardware execution
pass_manager.pre_init = ffsim.qiskit.PRE_INIT
isa_circuit = pass_manager.run(circuit)
print(f"Gate counts (w/ pre-init passes): {isa_circuit.count_ops()}")

Output:

Gate counts (w/o pre-init passes): OrderedDict({'rz': 7579, 'sx': 6106, 'ecr': 2316, 'x': 336, 'measure': 32, 'barrier': 1})
Gate counts (w/ pre-init passes): OrderedDict({'rz': 4088, 'sx': 3125, 'ecr': 1262, 'x': 201, 'measure': 32, 'barrier': 1})

3. Executar no hardware de destino

Após otimizar o circuito para execução em hardware, estamos prontos para executá-lo no hardware de destino e coletar amostras para a estimativa da energia do estado fundamental. Como temos apenas um circuito, vamos usar o modo de execução de tarefas qiskit-ibm-runtime e executar nosso circuito.

from qiskit_ibm_runtime import SamplerV2 as Sampler

sampler = Sampler(mode=backend)
sampler.options.dynamical_decoupling.enable = True

job = sampler.run([isa_circuit], shots=10_000)  # Takes approximately 5sec of QPU time

# Run cell after IQX job completion
primitive_result = job.result()
pub_result = primitive_result[0]
counts = pub_result.data.meas.get_counts()

4. Resultados pós-processamento

A parte de pós-processamento do fluxo de trabalho do SQD pode ser resumida no diagrama a seguir.

A amostragem do ansatz LUCJ na base computacional gera um conjunto de configurações com ruído $\tilde{\mathcal{\chi}}$, que são usadas na rotina de pós-processamento. Isso envolve um método chamado (detalhes discutidos posteriormente) de recuperação de configuração para corrigir probabilisticamente as configurações com números incorretos de elétrons. As configurações apenas com números de elétrons corretos $\tilde{\mathcal{\chi}}_{R}$ são então subamostradas e distribuídas em vários lotes com base na frequência de aparecimento de cada configuração exclusiva. Cada lote de amostras define um subespaço ( $\mathcal{S^{(k)}}$ ). Em seguida, o Hamiltoniano molecular, $H$, é projetado em subespaços:

Cada Hamiltoniano projetado $H_{\mathcal{S}^{(k)}}$ é então alimentado em um Eigensolver, onde é diagonalizado para calcular valores e vetores próprios para reconstruir um estado próprio. Nesta lição, projetamos e diagonalizamos o Hamiltoniano usando o pacote qiskit-addon-sqd , que utiliza o método de Davidson do site PySCF para a diagonalização.

Em seguida, coletamos o menor valor próprio (energia) dos lotes e também calculamos a ocupação orbital média, $\text{n}$. As informações de ocupação média são usadas na etapa de recuperação de configuração para corrigir probabilisticamente as configurações de ruído.

Em seguida, explicamos detalhadamente o loop de recuperação de configuração autoconsistente e mostramos exemplos concretos de código para implementar as etapas mencionadas acima para estimar a energia do estado fundamental do hamiltoniano $N_2$.

4.1 Recuperação da configuração: visão geral

Cada bit em uma cadeia de bits (determinante de Slater) representa um orbital de spin. A metade direita de uma cadeia de bits representa orbitais de spin para cima e a metade esquerda representa orbitais de spin para baixo. Um 1 significa que o orbital está ocupado por um elétron, e um 0 significa que o orbital está vazio. Sabemos a priori o número correto de partículas (elétron de spin ascendente e elétron de spin descendente). Suponha que tenhamos um determinante $x$ com $N_x$ elétrons (ou seja, há $N_x$ números de $1$ s na cadeia de bits) nele. O número correto de partículas é $N$. Se $N_x \neq N$, então sabemos que a cadeia de bits está corrompida por ruído. A rotina de configuração autoconsistente tenta corrigir a cadeia de bits invertendo probabilisticamente $|N_x - N|$ bits, aproveitando as informações de ocupação orbital média. A ocupação orbital média ( $n$ ) nos informa a probabilidade de um orbital ser ocupado por um elétron. Se $N_x < N$, temos menos elétrons e precisamos inverter alguns $0$ s para $1$ s e vice-versa.

A probabilidade de inversão pode ser $|x[i] - avg\_occupancy[i]|$ para i-th spin orbital. Em [2], os autores usaram uma probabilidade ponderada de inversão usando a função ReLU modificada.

Aqui, $h$ define o local do "canto" da função ReLU, e o parâmetro $\delta$ define o valor da função ReLU no canto. Para $\delta = 0$, $w$ se torna a função ReLU verdadeira, e para $\delta >0$, se torna ReLU* modificada*. No artigo, os autores usaram $\delta = 0.01$ e $h =$ número de partículas alfa (ou beta)/número de orbitais de spin alfa (ou beta) $= N/M$ (fator de preenchimento).

A ocupação orbital média ( $n$ ) não é conhecida a priori. A primeira iteração da estimativa do estado fundamental começa com configurações com apenas números de partículas corretos em ambas as espécies de spin. Após a primeira iteração, temos uma estimativa do estado fundamental e, usando a estimativa, podemos construir o primeiro palpite de $n$. Esse palpite de $n$ é usado para recuperar configurações, executar a próxima iteração de estimativa do estado fundamental e refinar de forma autoconsistente o palpite de $n$. O processo se repete até que um critério de parada seja atendido.

Considere o seguinte exemplo para $N = 2$ e $x = |1000\rangle$ ( $N_x = 1$ ). Precisamos inverter um dos 0s para 1 a fim de corrigi-lo para números de partículas, e as opções são 1100, 1010 e 1001. Com base na probabilidade de inversão, uma das opções será selecionada como configuração recuperada (ou a cadeia de bits com o número correto de partículas).

Suponha que, na primeira iteração, executemos dois lotes e que os estados de solo estimados a partir deles sejam:

Usando os estados da base computacional e suas amplitudes, podemos calcular a probabilidade de ocupações de elétrons (em resumo, ocupações ) por spin-orbital (qubit) (observe que a probabilidade = |amplitude| $^2$ ). Abaixo, tabulamos as ocupações por qubit para cada bitstring que aparece no estado fundamental estimado e calculamos a ocupação orbital total para um lote. Observe que, de acordo com a convenção de ordenação do Qiskit, o bit mais à direita representa qubit-0 ( Q0 ), e o bit mais à esquerda representa Q3.

Ocupação ( Batch0 ):

Ocupação ( Batch1 )

Ocupação (média dos lotes)

Usando a ocupação orbital média calculada acima, podemos encontrar as probabilidades de inversão para diferentes orbitais na configuração $x = \vert 1000 \rangle$. Como o orbital representado por Q3 já está ocupado e não precisa ser invertido, definimos seu p(flip) como $0$. Para os orbitais restantes, que estão desocupados, a probabilidade de inversão é $\vert x[i] - \text{n}[i] \vert$ cada. Junto com p(flip), também calculamos o peso da probabilidade associada ao flipping usando a função ReLU modificada descrita acima.

Probabilidade de flip ( $x = \vert 1000 \rangle$, $\delta = 0.01$, $h = N/M = 2/4 = 0.50$ )

Por fim, usando as probabilidades ponderadas acima, podemos inverter um dos orbitais desocupados Q2, Q1 e Q0. Com base nos valores acima, o Q0 provavelmente será invertido, e uma possível configuração recuperada pode ser $\vert \text{1001} \rangle$.

O processo completo de recuperação de configuração autoconsistente pode ser resumido da seguinte forma:

Primeira iteração: Suponha que as cadeias de bits (configurações ou determinantes de Slater) geradas pelo computador quântico formem um conjunto $\widetilde{\chi}$, que inclui configurações com o número correto ( $\widetilde{\chi}_{correct}$ ) e incorreto ( $\widetilde{\chi}_{incorrect}$ ) de partículas em cada setor de spin.

1. As configurações de ( $\widetilde{\chi}_{correct}$ ) são amostradas aleatoriamente para criar lotes $(\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)})$ de vetores para projeção de subespaço. O número de lotes e amostras em cada lote são parâmetros definidos pelo usuário. Quanto maior o número de amostras em cada lote, maior a dimensão do subespaço e mais exigente em termos computacionais se torna a diagonalização. Por outro lado, um número muito pequeno de amostras pode perder os vetores de suporte do estado básico e levar a uma estimativa incorreta.
2. Execute o solucionador de estado próprio (ou seja, projeção no subespaço e diagonalização) nos lotes e obtenha estados próprios aproximados. $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
3. A partir dos estados próprios aproximados, construa a primeira estimativa para $n$.

Iterações subsequentes:

1. Usando $n$, corrija as configurações com o número de partículas errado em $\widetilde{\chi}_{incorrect}$. Suponha que as chamemos de $\widetilde{\chi}_{correct\_new}$. Então, $\widetilde{\chi}_{recovered} (\widetilde{\chi}_{R}) = \widetilde{\chi}_{correct} \cup \widetilde{\chi}_{correct\_new}$ forma o novo conjunto de configurações com números de partículas corretos.
2. $\widetilde{\chi}_{R}$ é amostrado para criar lotes $\mathcal{S}^{(1)}, \cdots, \mathcal{S}^{(K)}$.
3. O solucionador de estado próprio é executado com novos lotes e gera novas estimativas de estados fundamentais $|\psi^{(1)}\rangle, \cdots, |\psi^{(K)}\rangle$.
4. A partir dos estados próprios aproximados, construa uma estimativa refinada para $n$.
5. Se o critério de parada não for atendido, volte para a etapa 2.1.

4.2 Estimativa do estado fundamental

Primeiro, transformaremos as contagens em uma matriz de bitstring e em uma matriz de probabilidade para pós-processamento.

Cada linha da matriz representa uma bitstring exclusiva. Como os qubits são indexados a partir da direita de uma cadeia de bits no Qiskit, a coluna 0 representa o qubit N-1 e a coluna N-1 representa o qubit 0, em que N é o número de qubits.

Os orbitais alfa são representados no intervalo de índice de coluna (N, N/2] (metade direita), e os orbitais beta são representados no intervalo de coluna (N/2, 0] (metade esquerda).

from qiskit_addon_sqd.counts import counts_to_arrays

# Convert counts into bitstring and probability arrays
bitstring_matrix_full, probs_arr_full = counts_to_arrays(counts)

Há algumas opções controladas pelo usuário que são importantes para essa técnica:

• iterations: Número de iterações de recuperação de configuração autoconsistente
• n_batches: Número de lotes de configurações usados pelas diferentes chamadas ao solucionador de estado próprio
• samples_per_batch: Número de configurações exclusivas a serem incluídas em cada lote
• max_davidson_cycles: Número máximo de ciclos Davidson executados por cada eigensolver

import numpy as np
from qiskit_addon_sqd.configuration_recovery import recover_configurations
from qiskit_addon_sqd.fermion import (
    bitstring_matrix_to_ci_strs,
    solve_fermion,
)
from qiskit_addon_sqd.subsampling import postselect_and_subsample

rng = np.random.default_rng(24)
# SQD options
iterations = 5

# Eigenstate solver options
n_batches = 5
samples_per_batch = 500
max_davidson_cycles = 300

# Self-consistent configuration recovery loop
e_hist = np.zeros((iterations, n_batches))  # energy history
s_hist = np.zeros((iterations, n_batches))  # spin history
occupancy_hist = []
avg_occupancy = None
for i in range(iterations):
    print(f"Starting configuration recovery iteration {i}")
    # On the first iteration, we have no orbital occupancy information from the
    # solver, so we begin with the full set of noisy configurations.
    if avg_occupancy is None:
        bs_mat_tmp = bitstring_matrix_full
        probs_arr_tmp = probs_arr_full

    # If we have average orbital occupancy information, we use it to refine
    # the full set of noisy configurations.
    else:
        bs_mat_tmp, probs_arr_tmp = recover_configurations(
            bitstring_matrix_full,
            probs_arr_full,
            avg_occupancy,
            num_elec_a,
            num_elec_b,
            rand_seed=rng,
        )

    # Create batches of subsamples. We postselect here to remove configurations
    # with incorrect hamming weight during iteration 0, since no config recovery was performed.
    batches = postselect_and_subsample(
        bs_mat_tmp,
        probs_arr_tmp,
        hamming_right=num_elec_a,
        hamming_left=num_elec_b,
        samples_per_batch=samples_per_batch,
        num_batches=n_batches,
        rand_seed=rng,
    )

    # Run eigenstate solvers in a loop. This loop should be parallelized for larger problems.
    e_tmp = np.zeros(n_batches)
    s_tmp = np.zeros(n_batches)
    occs_tmp = []
    coeffs = []
    for j in range(n_batches):
        strs_a, strs_b = bitstring_matrix_to_ci_strs(batches[j])
        print(f"  Batch {j} subspace dimension: {len(strs_a) * len(strs_b)}")
        energy_sci, coeffs_sci, avg_occs, spin = solve_fermion(
            batches[j],
            hcore,
            eri,
            open_shell=open_shell,
            spin_sq=spin_sq,
            max_cycle=max_davidson_cycles,
        )
        energy_sci += nuclear_repulsion_energy
        e_tmp[j] = energy_sci
        s_tmp[j] = spin
        occs_tmp.append(avg_occs)
        coeffs.append(coeffs_sci)

    # Combine batch results
    avg_occupancy = tuple(np.mean(occs_tmp, axis=0))

    # Track optimization history
    e_hist[i, :] = e_tmp
    s_hist[i, :] = s_tmp
    occupancy_hist.append(avg_occupancy)

Output:

Starting configuration recovery iteration 0
  Batch 0 subspace dimension: 21609
  Batch 1 subspace dimension: 21609
  Batch 2 subspace dimension: 21609
  Batch 3 subspace dimension: 21609
  Batch 4 subspace dimension: 21609
Starting configuration recovery iteration 1
  Batch 0 subspace dimension: 609961
  Batch 1 subspace dimension: 616225
  Batch 2 subspace dimension: 627264
  Batch 3 subspace dimension: 633616
  Batch 4 subspace dimension: 624100
Starting configuration recovery iteration 2
  Batch 0 subspace dimension: 564001
  Batch 1 subspace dimension: 605284
  Batch 2 subspace dimension: 582169
  Batch 3 subspace dimension: 559504
  Batch 4 subspace dimension: 591361
Starting configuration recovery iteration 3
  Batch 0 subspace dimension: 550564
  Batch 1 subspace dimension: 549081
  Batch 2 subspace dimension: 531441
  Batch 3 subspace dimension: 527076
  Batch 4 subspace dimension: 531441
Starting configuration recovery iteration 4
  Batch 0 subspace dimension: 544644
  Batch 1 subspace dimension: 580644
  Batch 2 subspace dimension: 527076
  Batch 3 subspace dimension: 531441
  Batch 4 subspace dimension: 537289

4.3 Discussão dos resultados

O primeiro gráfico mostra que, após algumas iterações, estimamos a energia do estado fundamental dentro de ~24 mH (a precisão química é normalmente aceita como sendo de 1 kcal/mol $\approx$ 1.6 mH ). O segundo gráfico mostra a ocupação média de cada orbital espacial após a iteração final. Podemos ver que tanto os elétrons de spin para cima quanto os de spin para baixo ocupam os primeiros cinco orbitais com alta probabilidade em nossas soluções.

Embora a energia estimada do estado fundamental seja decente, ela não está dentro do limite de precisão química ( $\pm \approx 1.6$ mH ). Essa lacuna pode ser atribuída à pequena dimensão do subespaço que usamos acima para projeção e diagonalização. Como usamos samples_per_batch=500, o subespaço é abrangido por no máximo $500$ vetores, que são os vetores ausentes do suporte do estado básico. O aumento do parâmetro samples_per_batch deve melhorar a precisão à custa de mais recursos de computação clássicos e tempo de execução.

# Data for energies plot
x1 = range(iterations)
min_e = [np.min(e) for e in e_hist]
e_diff = [abs(e - exact_energy) for e in min_e]
yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5]

# Chemical accuracy (+/- 1 milli-Hartree)
chem_accuracy = 0.001

# Data for avg spatial orbital occupancy
y2 = occupancy_hist[-1][0] + occupancy_hist[-1][1]
x2 = range(len(y2))

import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(12, 6))

# Plot energies
axs[0].plot(x1, e_diff, label="energy error", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(yt1)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-6)
axs[0].axhline(
    y=chem_accuracy, color="#BF5700", linestyle="--", label="chemical accuracy"
)
axs[0].set_title("Approximated Ground State Energy Error vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy Error (Ha)", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f} Ha")
print(f"SQD energy: {min_e[-1]:.5f} Ha")
print(f"Absolute error: {e_diff[-1]:.5f} Ha")
plt.tight_layout()
plt.show()

Output:

Exact energy: -109.04667 Ha
SQD energy: -109.02234 Ha
Absolute error: 0.02434 Ha

Exercício para o leitor

Aumente progressivamente o parâmetro samples_per_batch (por exemplo, de $1000$ para $10000$ em uma etapa de $1000$; permitido pela memória do seu computador) e compare as energias estimadas do estado fundamental.

Referências

[1] M. Motta et al, "Unindo intuição física e eficiência de hardware para estados eletrônicos correlacionados: o cluster unitário local Jastrow ansatz para estrutura eletrônica" (2023). Química. Ciências, 2023, 14, 11213.

[2] J. Robledo-Moreno et al, "Química além das soluções exatas em um supercomputador centrado em quantum" (2024). arXiv:quant-ph/2405.05068.