Adiabatic AF state preparation

[ComputationalBasis]

Hello everyone,

I am currently trying to simulate the adiabatic antiferromagnetic state preparation. Here is the code I made:

import torch
from torch import tensor
from qadence import sample, TimeParameter, FeatureParameter, X, Y, Z, N, add, kron, Register, HamEvo, run, BackendName, PI, expectation
from qadence.analog.constants import C6_DICT
from pyqtorch.utils import SolverType
import numpy as np
from sympy import cos, sin, exp, tanh

ode_solver = SolverType.KRYLOV_SE  # time-dependent Schrodinger equation solver method
n_steps_hevo = 10000  # integration time steps used by solver

# Parameters
c_6 = C6_DICT[60] # [rad*μm^6/μs]
U = 2*PI # [rad/μs]
omega_max = 2*U
delta_0 = -6*U
delta_f = 2*U

# Time parameters on pulses Omega and delta
t_rise = 0.252 # [μs]
t_fall = 0.5 # [μs]
t_total = 1.552 # [μs]

# Register
n = 3
dx = (c_6 / U) ** (1/6)  # [μm]
reg = Register.rectangular_lattice(qubits_row = 3, qubits_col = 3, spacing=dx)


# Initial state: ground state
initial_state = torch.zeros(2**(n**2), dtype=torch.complex128).unsqueeze(0)
initial_state[0] = 1

# Hamiltonian parameters

t = TimeParameter("t")
omega_param = FeatureParameter("omega") # [rad/μs]
delta_param = FeatureParameter("delta") # [rad/μs]

a = 50000
def f(t):
    return 0.5 * (tanh(a * t) + 1)

omega_t = omega_max/t_rise * t * (f(t) - f(t-t_rise)) + omega_max * (f(t-t_rise) - f(t-(t_total-t_fall))) + (omega_max/t_fall * (t_total-t)) * (f(t-(t_total-t_fall)) - f(t-t_total))
delta_t = delta_0 * (f(t+0.1) - f(t-t_rise)) + (((delta_f-delta_0)*t/(t_total-t_fall-t_rise)) + delta_0 - (delta_f-delta_0)*t_rise/(t_total-t_fall-t_rise)) * (f(t-t_rise) - f(t-(t_total-t_fall)))  + delta_f * (f(t-(t_total-t_fall)) - f(t-(t_total+0.1)))


Rb = (c_6/omega_t)**(1/6) 
cociente = Rb/dx

# Hamiltonian
h_x = (omega_param * omega_t / 2) * add(X(i) for i in range(n))
h_n = -1.0 * delta_param * delta_t * add(N(i) for i in range(n))

h_d = h_x + h_n  # Driving Hamiltonian

h_nn = U * (N(0)@N(1) + N(0)@N(1) + N(1)@N(2)
          + N(2)@N(5) + N(5)@N(8) + N(8)@N(7)
          + N(7)@N(6) + N(6)@N(3) + N(3)@N(0)
          + N(4)@N(1) + N(4)@N(3) + N(4)@N(7)
          + N(4)@N(5))# Ising

h = h_d + h_nn

# Duration time evolution
duration = t_total * 1e3  #  [ns] 1e3 ns = 1 μs


# Omega and delta values
param_values = {
    "omega": tensor(1),
    "delta": tensor(1)
}

# Evolution
evo = HamEvo(h, parameter=t, duration=duration/1000)

# Sample
xs = sample(reg, evo, state=initial_state, values=param_values, n_shots=10000)

The thing is, I used the exact same pulses Omega and delta that are being used on this pulser tutorial (you can check the pulses by running the omega_t and delta_t lines) and the same Rydberg coefficient. I also used the exact same register:

The issue is that I am not getting the AF state whatsoever. The result I obtained is:

Con 10000 shots: [OrderedCounter({'010100001': 85, '010100010': 81, '010100100': 77, '010100110': 74, '010100011': 69, '010100101': 65, '010100000': 64, '001011101': 62, '111101100': 62, '001010000': 61, '001011011': 61, '111101000': 60, '001011100': 59, '001011111': 59, '111000011': 59, '010100111': 58, '001011000': 57, '111000100': 56, '001011010': 55, '111000110': 55, '001010001': 53, '100110010': 53, '111000010': 53, '111000111': 53, '111101110': 53, '100110111': 52, '111000001': 52, '001010010': 51, '001010011': 51, '010001001': 51, '010110101': 51, '111101001': 51, '001011110': 50, '111000101': 50, '001010111': 49, '010001100': 49, '111000000': 49, '001011001': 48, '111111111': 48, '100010000': 47, '100010011': 47, '010001110': 46, '010001111': 46, '100010111': 46, '100110000': 45, '111101101': 45, '010001101': 44, '010110110': 44, '100110001': 44, '111101010': 44, '001000100': 43, '001010110': 43, '010001010': 43, '100010010': 43, '100110011': 43, '111111001': 43, '010001011': 42, '100010101': 42, '100010110': 42, '111111010': 42, '111111101': 42, '010001000': 41, '100010100': 41, '111101011': 41, '001000011': 40, '100110101': 40, '111011001': 40, '111110111': 40, '111111011': 40, '010110001': 39, '010110100': 39, '011111010': 39, '100110110': 39, '001000111': 38, '001010100': 38, '001010101': 38, '010110010': 38, '000111011': 37 ....

Maybe the problem lies when defining the Hamiltonian, as I do not exactly understand the Parameter thing. Can someone check if the Hamiltonian is well-defined? Also, maybe I have to drop this idea as the code takes too long to compute, if anyone has a better idea of study with Qadence I am open to hear it.

[ComputationalBasis]

Hello again everyone,

Somehow I managed to fix the code shown above (I'm impressed by myself), if anyone wants to see it, here it is:

import torch
from torch import tensor
from qadence import sample, product_state, TimeParameter, FeatureParameter, X, Y, Z, N, add, kron, Register, HamEvo, run, BackendName, PI, expectation
from qadence.analog.constants import C6_DICT
from pyqtorch.utils import SolverType
import numpy as np
from sympy import cos, sin, exp, tanh

ode_solver = SolverType.KRYLOV_SE  # time-dependent Schrodinger equation solver method
n_steps_hevo = 10000  # integration time steps used by solver

# Parameters
c_6 = C6_DICT[60] # [rad*μm^6/μs]
U = 2*PI # [rad/μs]

# Time parameters on pulses Omega and delta
t_rise = 0.252 # [μs]
t_fall = 0.5 # [μs]
t_total = 1.552 # [μs]

# Register
n = 3
dx = (c_6 / U) ** (1/6)  # [μm]
reg = Register.rectangular_lattice(qubits_row=n, qubits_col=n, spacing=dx)


# Initial state: ground state
initial_state = product_state("000000000")

# Hamiltonian parameters

t = TimeParameter("t")
omega_max_param = FeatureParameter("omega_max") # [rad/μs]
delta_0_param = FeatureParameter("delta_0") # [rad/μs]
delta_f_param = FeatureParameter("delta_f") # [rad/μs]

a = 50000
def f(t):
    return 0.5 * (tanh(a * t) + 1)

omega_t = omega_max_param/t_rise * t * (f(t) - f(t-t_rise)) + omega_max_param * (f(t-t_rise) - f(t-(t_total-t_fall))) + (omega_max_param/t_fall * (t_total-t)) * (f(t-(t_total-t_fall)) - f(t-t_total))
delta_t = delta_0_param * (f(t+0.1) - f(t-t_rise)) + (((delta_f_param-delta_0_param)*t/(t_total-t_fall-t_rise)) + delta_0_param - (delta_f_param-delta_0_param)*t_rise/(t_total-t_fall-t_rise)) * (f(t-t_rise) - f(t-(t_total-t_fall)))  + delta_f_param * (f(t-(t_total-t_fall)) - f(t-(t_total+0.1)))


# Hamiltonian
h_x = (omega_t / 2) * add(X(i) for i in range(3*n))
h_n = -1.0 * delta_t * add(N(i) for i in range(3*n))

h_d = h_x + h_n  # Driving Hamiltonian
h_nn = U * (N(0)@N(1) + N(1)@N(2)
          + N(2)@N(5) + N(5)@N(8) + N(8)@N(7)
          + N(7)@N(6) + N(6)@N(3) + N(3)@N(0)
          + N(4)@N(1) + N(4)@N(3) + N(4)@N(7)
          + N(4)@N(5)) # Ising

h = h_d + h_nn

# Duration time evolution
duration = t_total * 1e3  #  [ns] 1e3 ns = 1 μs


# Omega and delta values
param_values = {
    "omega_max": tensor(2*U),
    "delta_0": tensor(-6*U),
    "delta_f": tensor(2*U)
}

# Evolution
evo = HamEvo(h, parameter=t, duration=duration/1000)

# Sample
xs = sample(reg, evo, state=initial_state, values=param_values, n_shots=10000)

which gives the desired AF state:

Con 10000 shots: [OrderedCounter({'101010101': 2569, '101011101': 336, '101110101': 336, '101010111': 317, '101100011': 217, '101001101': 204, '101000101': 198, '101001110': 191, '101000111': 183, '101100101': 180, '001110101': 132, '101100111': 126, '101101010': 124, '101101101': 116, '101001111': 111, '011100101': 105, '101010100': 104, '101101011': 104, '101110001': 101, '001010101': 99, '100011101': 98, '101101110': 94, '010101101': 89, '101110011': 88, '101010001': 86, '101010011': 85, '101011100': 83, '101010110': 80, '011100011': 79, '101101111': 76, '101011110': 69, '100010101': 67, '101011111': 65, '011100111': 64, '001100101': 63, '011010101': 63, '101110111': 61, '011101101': 56, '101100001': 56, '100001101': 55, '011110101': 54, '101000011': 52, '001100011': 51, '100101101': 50 ...

@jpmoutinho you can check if this is good enough, I am happy with the result :)

[ComputationalBasis]

Additionaly, if anyone is curious, this is the evolution of the whole system under the hamiltonian shown in the selected reply:

the phase transition from the ground state to the antiferromagnetic state is clearly visible, although the qubit 4 is not that populated since it is the one in the middle of the square lattice, and therefore it is interacting with 4 nearest neighbours.

EDIT: there was a slight error in the code. In the interacting Hamiltonian, the term N(0)@N(1) was repeated twice. That is why there is a slight asymmetry with the figure shown above. If corrected, all 4 atoms at the corners have identical evolution over time.