Gere circuitos de evolução temporal de Krylov (SKQD)

Importante

Os conceitos deste guia estão atualmente disponíveis apenas na API Python. Funcionalidade equivalente será disponibilizada na API C em uma versão futura.

A diagonalização quântica de Krylov baseada em amostras (SKQD) é uma variante quântico-cêntrica da diagonalização de subespaço. Em vez de amostrar um único ansatz variacional (como na diagonalização quântica baseada em amostras simples (SQD)) ou um conjunto de circuitos de evolução temporal aleatorizados (como no SqDRIFT), ela amostra uma base de Krylov: a família de estados

$$\lvert \psi_k \rangle = e^{-i k \, \Delta t \, H} \lvert \psi_0 \rangle, \qquad k = 0, 1, \dots, D - 1,$$

obtidos pela evolução de um estado de referência $\lvert \psi_0 \rangle$ por múltiplos crescentes de um passo de tempo fixo $\Delta t$. As bitstrings amostradas desses $D$ circuitos povoam o subespaço de Krylov no qual o Hamiltoniano é, em seguida, diagonalizado classicamente.

Este guia reproduz a etapa de construção de circuitos do algoritmo SKQD para o modelo de Anderson de impureza única (SIAM), uma impureza magnética acoplada a um banho sem interação, seguindo a publicação do SKQD.

Nota

A etapa de amostragem do SKQD só funciona se o estado fundamental for esparso na base computacional: a diagonalização clássica precisa ser capaz de reconstruí-lo a partir de um número gerenciável de bitstrings amostradas. Como mostrado a seguir, essa esparsidade não é automática; ela é uma propriedade da base de partícula única na qual o problema é expresso. Acertar a base é o ponto central desta construção.

1. Construa o Hamiltoniano do SIAM como um operador fermiônico

O modelo de Anderson de impureza única coloca um orbital de impureza com interação (índice $0$) em contato com um banho sem interação de $L$ sítios dispostos em cadeia. Seu Hamiltoniano se divide em uma parte de um corpo $H_1$ e uma parte de dois corpos $H_2$,

$$H = H_1 + H_2,$$

com a parte de um corpo reunindo a energia local da impureza (o potencial químico $\mu$), a hibridização impureza–banho $V$ e o salto entre vizinhos mais próximos do banho $t$,

$$H_1 = \mu \sum_\sigma a^\dagger_{0\sigma} a_{0\sigma} - V \sum_\sigma \left( a^\dagger_{0\sigma} a_{1\sigma} + \text{h.c.} \right) - t \sum_{\sigma} \sum_{i=1}^{L-1} \left( a^\dagger_{i\sigma} a_{i+1,\sigma} + \text{h.c.} \right),$$

e a parte de dois corpos carregando a repulsão de Coulomb local $U$ apenas na impureza,

$$H_2 = U \, a^\dagger_{0\uparrow} a_{0\uparrow} a^\dagger_{0\downarrow} a_{0\downarrow}.$$

Este guia trabalha no ponto de simetria partícula-buraco $\mu = -U/2$. Essas duas partes correspondem diretamente às integrais de um corpo h1e (salto do banho, hibridização da impureza, energia local da impureza) e à única integral de dois corpos h2e (o termo de Coulomb da impureza). Primeiro, monte-as na base de posição, na qual o modelo é naturalmente definido.

>>> import numpy as np
>>>
>>> norb = 8
>>> nelec = (norb // 2, norb // 2)  # half filling
>>> hopping, onsite, hybridization = 1.0, 10.0, 1.0
>>> chemical_potential = -0.5 * onsite
>>>
>>> # one-body integrals: bath hopping, impurity hybridization, and impurity on-site energy
>>> h1e = np.zeros((norb, norb))
>>> np.fill_diagonal(h1e[:, 1:], -hopping)
>>> np.fill_diagonal(h1e[1:, :], -hopping)
>>> h1e[0, 1] = h1e[1, 0] = -hybridization
>>> h1e[0, 0] = chemical_potential
>>>
>>> # two-body integrals: on-site Coulomb repulsion on the impurity orbital
>>> h2e = np.zeros((norb,) * 4)
>>> h2e[0, 0, 0, 0] = onsite

2. Mude para a base de momento

A parte do banho no SIAM é uma cadeia de férmions livres, de modo que é diagonalizada por uma mudança para a base de momento. Esta é a etapa decisiva para o SKQD. Como o banho é invariante por translação, o estado fundamental não é esparso na base de posição, mas é esparso na base de momento, na qual o Hamiltoniano do banho é diagonal. A rotação de orbitais de partícula única que realiza essa mudança é álgebra linear específica do problema. Ela diagonaliza a matriz de salto do banho sem alterar o orbital da impureza e, em seguida, realoca a impureza para um sítio central:

>>> def momentum_basis(norb):
...     """Orbital rotation from the position to the momentum basis."""
...     n_bath = norb - 1
...     hopping_matrix = np.zeros((n_bath, n_bath))
...     np.fill_diagonal(hopping_matrix[:, 1:], -1)
...     np.fill_diagonal(hopping_matrix[1:, :], -1)
...     _, vecs = np.linalg.eigh(hopping_matrix)
...     orbital_rotation = np.zeros((norb, norb))
...     orbital_rotation[0, 0] = 1.0
...     orbital_rotation[1:, 1:] = vecs
...     new_index = n_bath // 2
...     perm = np.r_[1:(new_index + 1), 0, (new_index + 1):norb]
...     return orbital_rotation[:, perm]
...
>>> orbital_rotation = momentum_basis(norb)

Ao contrário de um fluxo de trabalho variacional, aqui as integrais do Hamiltoniano são rotacionadas para a base de momento, de modo que o circuito inteiro (Hamiltoniano e estado) fica na base em que o estado fundamental é esparso e, portanto, as bitstrings amostradas são informativas. Rotacionar as integrais de um e de dois corpos é uma mudança de base padrão dos tensores de estrutura eletrônica:

>>> rotation = orbital_rotation.T.conj()
>>> h1e_momentum = np.einsum("ab,Aa,Bb->AB", h1e, rotation, rotation.conj(), optimize="greedy")
>>> h2e_momentum = np.einsum(
...     "abcd,Aa,Bb,Cc,Dd->ABCD",
...     h2e,
...     rotation,
...     rotation.conj(),
...     rotation,
...     rotation.conj(),
...     optimize="greedy",
... )

Como a rotação mantém o orbital da impureza fixo (ela atua apenas no banho), a interação local permanece inalterada. h2e_momentum continua sendo um único termo número-número, assim como h2e. A única mudança é que a impureza foi realocada para um índice central. Isso é visível diretamente nos tensores de integrais, que têm cada um uma entrada diferente de zero:

>>> print(f"h2e nonzero at:          {np.argwhere(np.abs(h2e) > 1e-12).tolist()}")
h2e nonzero at:          [[0, 0, 0, 0]]
>>> print(f"h2e_momentum nonzero at: {np.argwhere(np.abs(h2e_momentum) > 1e-12).tolist()}")
h2e_momentum nonzero at: [[3, 3, 3, 3]]

Essa localidade é o que mantém a interação barata de sintetizar mais adiante.

Monte as integrais de um corpo em um FermionOperator com o construtor de integrais eletrônicas from_1body_tril_spin_sym(), que espera que elas estejam empacotadas na ordenação triangular inferior. A parte de dois corpos é um único termo número-número; então, em vez de passá-la pela maquinaria de integrais de dois corpos, construa-a diretamente com cre()/ann(): o operador de Coulomb local $U \, n_{p\uparrow} n_{p\downarrow}$ para o modo de impureza $p$ (seu modo $\uparrow$ no índice p, seu modo $\downarrow$ no índice p + norb no layout de spin em blocos). Ambos os construtores duplicam o spin dos norb orbitais espaciais em 2 * norb modos fermiônicos.

Mantenha a parte de um corpo $H_1$ e a parte de dois corpos $H_2$ como operadores separados, em vez de somá-los imediatamente, de acordo com os termos $H_1$ e $H_2$ das equações da seção 1. Construir cada parte uma única vez (nas duas bases) permite reutilizá-las adiante. Somadas, elas formam o Hamiltoniano completo para a diagonalização exata em qualquer uma das bases. Individualmente, $H_2$ (base de momento) é o operador que o passo de Trotter evolui na etapa 4.

>>> from qiskit_fermions.operators import FermionOperator, ann, cre
>>>
>>> # U n_imp_up n_imp_down: the impurity's up/down modes sit at `imp` and `imp + norb` in the
>>> # block-spin layout (the impurity is orbital 0 in the position basis, orbital 3 in the
>>> # momentum basis -- the single relocated index found above)
>>> imp_pos, imp_mom = 0, 3
>>>
>>> one_body_position = FermionOperator.from_1body_tril_spin_sym(
...     h1e[np.tril_indices(norb)], norb
... )
>>> two_body_position = FermionOperator.from_dict(
...     {(cre(imp_pos), ann(imp_pos), cre(imp_pos + norb), ann(imp_pos + norb)): onsite}
... )
>>> one_body_momentum = FermionOperator.from_1body_tril_spin_sym(
...     h1e_momentum[np.tril_indices(norb)], norb
... )
>>> two_body_momentum = FermionOperator.from_dict(
...     {(cre(imp_mom), ann(imp_mom), cre(imp_mom + norb), ann(imp_mom + norb)): onsite}
... )
>>>
>>> hamiltonian_position = (
...     (one_body_position + two_body_position).normal_ordered().simplify(atol=1e-14)
... )
>>> hamiltonian_momentum = (
...     (one_body_momentum + two_body_momentum).normal_ordered().simplify(atol=1e-14)
... )

A mudança de base é uma transformação de similaridade unitária, portanto não altera o espectro. Confirme que o operador está correto comparando seu menor autovalor no setor $(\text{norb}, \text{nelec})$ com uma diagonalização exata. O FermionOperator expõe um LinearOperator do SciPy (apoiado por um kernel nativo de produto matriz-vetor de interação de configurações completa (FCI)), de modo que ele pode ser passado diretamente para scipy.sparse.linalg.eigsh().

>>> import ffsim
>>> import scipy.sparse.linalg
>>>
>>> linop = ffsim.linear_operator(hamiltonian_momentum, norb, nelec)
>>> reference_energy, ground_state = scipy.sparse.linalg.eigsh(linop, k=1, which="SA")
>>> reference_energy = reference_energy[0]
>>> ground_state = ground_state[:, 0]
>>> print(f"exact ground-state energy: {reference_energy:.6f}")
exact ground-state energy: -13.422492

O ganho da mudança de base é a esparsidade. Para torná-lo concreto, conte quantos determinantes da base computacional são necessários para capturar 99% do peso do estado fundamental e compare a base de momento com a base de posição. A energia é independente da base, mas o número de determinantes não é:

>>> def determinants_for(vec, fraction=0.99):
...     weights = np.sort(np.abs(vec) ** 2)[::-1]
...     return int(np.searchsorted(np.cumsum(weights), fraction) + 1)
...
>>> # the same Hamiltonian in the position basis has an identical spectrum, but a dense ground state
>>> position_linop = ffsim.linear_operator(hamiltonian_position, norb, nelec)
>>> position_energy, position_ground = scipy.sparse.linalg.eigsh(position_linop, k=1, which="SA")
>>>
>>> print(f"position-basis energy: {position_energy[0]:.6f}")
position-basis energy: -13.422492
>>> print(f"determinants for 99% weight: {determinants_for(position_ground[:, 0])} (position)")
determinants for 99% weight: 2563 (position)
>>> print(f"determinants for 99% weight: {determinants_for(ground_state)} (momentum)")
determinants for 99% weight: 21 (momentum)

Os dois operadores compartilham a mesma energia, mas o estado fundamental na base de posição se espalha por milhares de determinantes, enquanto a base de momento o concentra em cerca de 20, e é por isso que a amostragem do SKQD funciona na base de momento e não na base de posição.

3. Prepare o estado de referência

A referência do SKQD $\lvert \psi_0 \rangle$ é a superposição de todas as excitações dos elétrons mais próximos do nível de Fermi para os modos de momento vazios mais próximos. A publicação do SKQD a constrói com uma rede escrita à mão de XXPlusYYGates. Ela parece fortemente correlacionada (tem centenas de amplitudes não nulas na base computacional), mas uma rede de rotações de dois modos que conservam o número é uma operação de partícula única (gaussiana fermiônica), de modo que o estado que ela produz é um único determinante de Slater. Isso significa que você pode prepará-lo com uma única porta PrepareSlaterDeterminant, dada a rotação de partícula única que a rede implementa.

Como cada XXPlusYYGate atua como uma rotação de Givens em um par de modos, essa rotação pode ser construída diretamente na representação de partícula única (uma matriz de Givens por porta), em vez de emitir manualmente as portas de dois qubits. As portas formam um padrão de tijolos (brickwork) de rotações entre vizinhos mais próximos que se expande a partir do nível de Fermi em camadas crescentes, espalhando cada elétron próximo ao nível de Fermi pelos modos vazios logo acima dele:

>>> def fermi_level_rotation(norb, nocc, n_excited=3):
...     """Orbital rotation exciting the near-Fermi electrons into the nearest empty modes."""
...     cos, sin = np.cos(np.pi / 4), np.sin(np.pi / 4)
...     fermi = nocc - 1  # highest occupied mode; the Fermi level sits between it and mode `nocc`
...     # Nearest-neighbor Givens rotations fanning out from the Fermi level in expanding layers:
...     # layer `depth` couples the mode `depth` steps below the level to the one above it, mixing
...     # the `n_excited` highest-occupied modes with the empty modes just above them. A final
...     # rotation carries amplitude up into the topmost empty mode of the window.
...     network = [
...         (lower, lower + 1)
...         for depth in range(n_excited)
...         for lower in range(fermi - depth, fermi + depth + 1, 2)
...     ]
...     network.append((fermi + n_excited, fermi + n_excited + 1))
...     unitary = np.eye(norb, dtype=complex)
...     for i, j in network:
...         givens = np.eye(norb, dtype=complex)
...         givens[i, i] = givens[j, j] = cos
...         givens[i, j], givens[j, i] = sin, -sin
...         unitary = givens @ unitary
...     return unitary
...
>>> reference_rotation = fermi_level_rotation(norb, nelec[0])

Nota

A rotação acima é expressa na base de momento, correspondendo ao Hamiltoniano. Os elétrons próximos ao nível de Fermi são os modos de momento ocupados de maior energia, e a rede os excita para os modos vazios de menor energia. É por isso que a referência deve ser construída depois da mudança de base da etapa anterior, e não antes.

A preparação de determinantes de Slater é feita por setor de spin, já que cada setor tem seu próprio formato (elétrons, orbitais). Uma única porta PrepareSlaterDeterminant (uma ocupação, com os nocc modos mais baixos preenchidos, junto com a rotação de partícula única) prepara um setor; a mesma porta é aplicada uma vez nos modos alfa 0..norb e uma vez nos modos beta norb..2*norb, já que ambos os setores de spin compartilham aqui a mesma ocupação de referência e a mesma rotação:

>>> from qiskit_fermions.circuit.library import PrepareSlaterDeterminant
>>>
>>> occupation = [i < nelec[0] for i in range(norb)]
>>> reference_preparation = PrepareSlaterDeterminant(occupation, reference_rotation)

Essa porta declarativa substitui a rede de XXPlusYYGate codificada à mão; o transpilador gera a síntese das rotações de Givens.

4. Monte os circuitos de Krylov

Cada circuito de Krylov é a preparação do estado de referência seguida da evolução sob $e^{-i k \, \Delta t \, H}$. Divida o Hamiltoniano em suas partes de um corpo e de dois corpos, $H = H_1 + H_2$, e aplique Trotter à evolução sobre essa divisão. Essa é a escolha decisiva para a profundidade do circuito. A evolução de um corpo $e^{-i t H_1}$ é uma operação de férmions livres (gaussiana fermiônica), portanto é uma rotação de base de partícula única (uma OrbitalRotation com unitária $e^{-i t \, \mathbf{h}_1}$) e não algo que precise passar por Trotter de Pauli termo a termo. A evolução de dois corpos $e^{-i t H_2}$ atua apenas sobre o único termo de interação local. Construir o passo dessa forma permite que o transpilador emita um brickwork raso de rotações de Givens mais uma única porta de interação de dois qubits, a estrutura rasa que a publicação do SKQD constrói à mão.

Já temos a parte de dois corpos. Ela é o operador two_body_momentum construído na etapa 2 (o único termo local). A evolução de um corpo $e^{-i t H_1}$ não depende do spin, portanto é uma única OrbitalRotation por setor, com unitária $e^{-i t \, \mathbf{h}_1}$ aplicada uma vez a cada setor de spin, com a mesma estrutura paralela da preparação do estado de referência.

Cada circuito de Krylov evolui por um tempo total $k \, \Delta t$. Isso é realizado com um produto de Trotter de segunda ordem de $k$ passos de tamanho $\Delta t$, de modo que o erro por passo permanece fixo à medida que a dimensão de Krylov cresce. Cada passo coloca uma rotação completa de um corpo (uma por setor de spin) entre dois meios passos da evolução de dois corpos (barata, de termo único). A dimensão de Krylov $D$ é o número desses circuitos, com dim variando de 0 (apenas o estado de referência) até $D - 1$.

>>> import scipy.linalg
>>>
>>> from qiskit_fermions.circuit import FermionicCircuit
>>> from qiskit_fermions.circuit.library import Evolution, OrbitalRotation
>>>
>>> num_modes = 2 * norb
>>> time_step = 0.2
>>>
>>> half_exp_h2 = Evolution(num_modes, two_body_momentum, time_step / 2)
>>> full_exp_h1 = OrbitalRotation(scipy.linalg.expm(-1j * time_step * h1e_momentum))
>>>
>>> def krylov_circuit(dim):
...     # closes over the fermionic gates defined in the outer scope:
...     # `reference_preparation`, `half_exp_h2` and `full_exp_h1`
...     circuit = FermionicCircuit(num_modes)
...     circuit.append(reference_preparation, range(norb))
...     circuit.append(reference_preparation, range(norb, num_modes))
...     for _ in range(dim):  # second-order Trotter product of `dim` steps
...         circuit.append(half_exp_h2, circuit.modes)
...         circuit.append(full_exp_h1, range(norb))
...         circuit.append(full_exp_h1, range(norb, num_modes))
...         circuit.append(half_exp_h2, circuit.modes)
...     return circuit

Nota

Os circuitos construídos aqui ainda não têm medições, porque os mesmos circuitos fermiônicos também alimentam a simulação exata na etapa 6, e uma medição não é uma operação unitária. A etapa 5 as adiciona com FermionicCircuit.measure_all(), que lê o modo i no bit clássico i.

Prefira medir o circuito fermiônico em vez do transpilado. Uma medição fermiônica nomeia o modo que lê, de modo que o layout e o roteamento a levam junto com o qubit em que esse modo foi parar, e a bitstring amostrada é indexada por modo. Medir o circuito transpilado, em vez disso, deixa para você a tarefa de desfazer final_index_layout(), e retorna silenciosamente resultados permutados se você esquecer. Isso depende da codificação de Jordan-Wigner na base de ocupação usada pelo preset, em que cada modo é carregado por um qubit.

5. Transpile para circuitos de qubits

Transpilar um FermionicCircuit mapeia suas portas fermiônicas para portas de qubits. Isso usa o mapeamento de Jordan-Wigner jordan_wigner(). As portas OrbitalRotation são sintetizadas em um brickwork de rotações de Givens, e a Evolution de termo único é reduzida a uma única porta de interação de dois qubits.

Como toda porta fermiônica aqui tem uma síntese padrão, você pode usar diretamente o generate_preset_jw_pass_manager() pronto, em vez de montar as etapas manualmente. Quaisquer argumentos nomeados são repassados à etapa de qubits; passe optimization_level=0 para ter uma imagem fiel, sem otimização, da profundidade sintetizada. Gerar a família completa de circuitos de Krylov é então um loop sobre a dimensão de Krylov. Mantenha os FermionicCircuits não transpilados junto com os circuitos de qubits transpilados. Os fermiônicos alimentam a simulação exata (vetor de estado) na etapa 6, portanto devem permanecer sem medições; measure_all(inplace=False) retorna uma cópia com medições para transpilar e deixa o original intacto, e essa cópia é o que um backend executaria.

>>> from qiskit_fermions.transpiler.presets import generate_preset_jw_pass_manager
>>>
>>> pass_manager = generate_preset_jw_pass_manager(optimization_level=0)
>>>
>>> def krylov_circuits(krylov_dim):
...     fermionic_circuits, circuits = [], []
...     for dim in range(krylov_dim):
...         fermionic_circuit = krylov_circuit(dim)
...         fermionic_circuits.append(fermionic_circuit)
...         measured = fermionic_circuit.measure_all(inplace=False)
...         circuits.append(pass_manager.run(measured))
...     return fermionic_circuits, circuits
...
>>> krylov_dim = 5
>>> fermionic_circuits, circuits = krylov_circuits(krylov_dim)

A dimensão de Krylov aparece diretamente na profundidade de dois qubits dos circuitos. A preparação do estado de referência é um custo fixo, e cada potência de Krylov adicional acrescenta mais um passo de Trotter de segunda ordem do operador de evolução temporal.

>>> two_qubit = lambda instr: instr.operation.num_qubits == 2
>>> [circuit.depth(two_qubit) for circuit in circuits]
[4, 12, 21, 30, 39]

Nota

Dividir o Hamiltoniano mantém essas profundidades modestas. A evolução de um corpo é uma verdadeira rotação de orbitais, então é sintetizada em um brickwork de XXPlusYYGate em qubits adjacentes. Cada meio passo da interação local é uma única RZZGate que acopla os dois modos de spin da impureza; o termo barato que o produto de segunda ordem coloca nos dois lados da rotação. Entregar o Hamiltoniano inteiro a uma única porta Evolution faria, em vez disso, o Trotter de Pauli da parte de um corpo termo a termo, propagando longas cadeias de Z de Jordan-Wigner pelos acoplamentos impureza-banho (de longo alcance) e inflando a profundidade em uma ordem de grandeza.

O primeiro circuito prepara apenas o determinante de referência, enquanto o último carrega o maior número de passos de Trotter (quatro, neste caso):

>>> circuits[0].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

>>> circuits[-1].draw("mpl", fold=-1, measure_arrows=False)
<Figure size ... with 1 Axes>

6. Amostre bitstrings dos circuitos de Krylov

Em hardware, você executaria os circuitos transpilados e mediria. Esta é a versão sem ruído (vetor de estado) dessa etapa, mostrando as bitstrings que o SKQD coletaria. Simule cada FermionicCircuit não transpilado com ffsim.apply_unitary() e amostre o vetor de estado resultante com ffsim.sample_state_vector().

Uma porta PrepareSlaterDeterminant faz validar e depois rotacionar na simulação (ela verifica se o estado de entrada ocupa seu determinante de referência e, em seguida, aplica a rotação). Portanto, a entrada é o determinante de ocupação simples: os nocc modos mais baixos preenchidos por spin, sem rotação. A própria porta aplica reference_rotation.

>>> from collections import Counter
>>>
>>> reference = ffsim.slater_determinant(
...     norb, (list(range(nelec[0])), list(range(nelec[1])))
... )
>>>
>>> shots = 1000
>>> counts = Counter()
>>> for dim, fermionic_circuit in enumerate(fermionic_circuits):
...     statevector = ffsim.apply_unitary(reference, fermionic_circuit, norb=norb, nelec=nelec)
...     samples = ffsim.sample_state_vector(
...         statevector, norb=norb, nelec=nelec, shots=shots, seed=dim
...     )
...     counts += Counter(samples)

Reunindo as amostras de todos os cinco circuitos de Krylov, o suporte é minúsculo: algumas centenas de bitstrings distintas entre os $\binom{8}{4}^2 = 4900$ determinantes do setor (4, 4). Essa concentração, consequência direta da esparsidade na base de momento estabelecida na etapa 2, torna tratável a diagonalização clássica seguinte.

Plotar as contagens com qiskit.visualization.plot_histogram() torna a esparsidade visível. Algumas configurações dominam, lideradas pelo determinante de referência.

>>> from qiskit.visualization import plot_histogram
>>>
>>> plot_histogram(dict(counts), number_to_keep=25)
<Figure size ... with 1 Axes>

7. Diagonalize no subespaço amostrado

As bitstrings amostradas são a entrada da metade clássica do SKQD. O Hamiltoniano é projetado no subespaço gerado pelas configurações amostradas e diagonalizado nele. Isso é feito com o pacote qiskit-addon-sqd, cuja função diagonalize_fermionic_hamiltonian() executa o laço completo de diagonalização quântica baseada em amostras. Ela constrói um subespaço a partir de lotes das configurações amostradas, diagonaliza o Hamiltoniano nele e refina o subespaço iterativamente usando a recuperação de configurações: inverter ocupações para reparar configurações que violam a simetria conhecida do número de partículas. A recuperação de configurações foi projetada para desfazer o dano do ruído do hardware, que corrompe as bitstrings amostradas para o setor errado de número de partículas; a amostragem aqui é livre de ruído, então não há nada a reparar nesse sentido.

O solver consome as contagens amostradas como um BitArray, construído diretamente a partir das counts agrupadas do passo 6. Passe num_bits explicitamente: sem ele, a largura é inferida da bitstring amostrada mais larga, que é mais estreita que 2 * norb sempre que nenhuma amostra ocupa o último modo, e então o solver rejeita o array. O solver diagonaliza na mesma base de momento que o restante da construção usa, então ele recebe diretamente as integrais h1e_momentum e h2e_momentum do passo 2, sem reempacotamento. A energia recuperada coincide com a referência exata com precisão bem abaixo de um mili-Hartree, a partir de apenas algumas centenas de configurações amostradas. Como a energia do SQD é um limite superior variacional do estado fundamental verdadeiro, ela sempre fica igual ou acima de reference_energy (ainda disponível do passo 2):

>>> from qiskit.primitives import BitArray
>>> from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian
>>>
>>> bit_array = BitArray.from_counts(dict(counts), num_bits=2 * norb)
>>>
>>> # record each iteration's batch results so we can watch the energy converge
>>> result_history = []
>>> result = diagonalize_fermionic_hamiltonian(
...     h1e_momentum,
...     h2e_momentum,
...     bit_array,
...     samples_per_batch=100,
...     norb=norb,
...     nelec=nelec,
...     num_batches=3,
...     max_iterations=3,
...     symmetrize_spin=True,
...     callback=result_history.append,
...     seed=np.random.default_rng(24),
... )
>>> print(f"SQD estimate: {result.energy:.4f} (exact: {reference_energy:.4f})")
SQD estimate: -13.4224 (exact: -13.4225)

Visualize a execução: a convergência da energia ao longo das iterações e a ocupação média de cada orbital espacial no estado fundamental recuperado. A energia por iteração é a menor entre os lotes daquela iteração, e a ocupação soma as orbital_occupancies sobre os dois setores de spin. Mesmo sem ruído a corrigir, a energia diminui de iteração em iteração. Cada iteração diagonaliza dentro de num_batches lotes de apenas samples_per_batch configurações subamostradas do conjunto, e o refinamento iterativo devido ao carregamento de bitstrings melhora continuamente quais configurações caem nesses lotes, de modo que a energia recuperada melhora:

>>> import matplotlib.pyplot as plt
>>> from matplotlib.ticker import MaxNLocator
>>>
>>> min_energies = [min(batch, key=lambda r: r.energy).energy for batch in result_history]
>>> occupancies = np.sum(result.orbital_occupancies, axis=0)
>>>
>>> fig, axs = plt.subplots(1, 2, figsize=(12, 5))
>>> _ = axs[0].plot(range(len(min_energies)), min_energies, marker="o", label="SQD energy")
>>> _ = axs[0].axhline(reference_energy, color="gray", linestyle="--", label="exact energy")
>>> _ = axs[0].set_xlabel("iteration")
>>> _ = axs[0].set_ylabel("energy")
>>> _ = axs[0].set_title("Energy vs SQD iteration")
>>> _ = axs[0].xaxis.set_major_locator(MaxNLocator(integer=True))
>>> _ = axs[0].legend()
>>> _ = axs[1].bar(range(norb), occupancies)
>>> _ = axs[1].set_xlabel("spatial orbital")
>>> _ = axs[1].set_ylabel("average occupancy")
>>> _ = axs[1].set_title("Occupancy per spatial orbital")
>>> fig
<Figure size ... with 2 Axes>

Isso fecha o laço do SKQD. Os circuitos na base de momento dos passos 1–5 produzem uma amostra esparsa (passo 6), e a diagonalização baseada em amostras acima transforma essa amostra de volta em uma energia do estado fundamental.

Próximos passos

Este guia executou todo o pipeline do SKQD em um simulador de vetor de estado sem ruído. Em hardware, os circuits transpilados seriam executados e medidos no lugar da amostragem do vetor de estado no passo 6, e as contagens resultantes alimentariam a diagonalização do passo 7 sem alterações. Consulte a documentação do Qiskit para obter ajuda na execução de circuitos e os tutoriais do addon SQD para informações sobre o pós-processamento de diagonalização de subespaço, incluindo seu comportamento em amostras ruidosas, onde a recuperação de configurações faz mais trabalho.

Leia o guia de relação com o ffsim para entender por que ffsim.apply_unitary() e ffsim.sample_state_vector() funcionam nativamente nos circuitos fermiônicos deste pacote, e por que o setor de número de partículas fixo usado no passo 6 torna eficiente em amostras a esparsidade da base de momento estabelecida no passo 2.