Hi @CatalinaAlbornoz - here is a mardkown export of my run without the output cells to stick under word count limit for posts on the forum
import pennylane as qml
from jax import numpy as jnp
import jax
jax.config.update('jax_enable_x64', True)
jax.config.update("jax_platform_name", "cpu")
Tapering the molecular Hamiltonian
#symbols = ['He','H']
#geometry = jnp.array([[0.0000, 0.0000, 0.0],
# [0.0000, 0.0000, 1.40]])
#n_electrons = 2
#charge = 1
symbols = ['Li','H']
geometry = jnp.array([[0.0000, 0.0000, -1.5065],
[0.0000, 0.0000, 1.5065]])
n_electrons = 4
charge = 0
molecule = qml.qchem.Molecule(symbols, geometry, charge=charge)
H, qubits = qml.qchem.molecular_hamiltonian(molecule, mapping='jordan_wigner')#, active_electrons=2)# active_electrons=2, active_orbitals=3, mapping='parity')
H
print (qubits, len(H))
generators = qml.symmetry_generators(H)
paulixops = qml.paulix_ops(generators, qubits)
for idx, generator in enumerate(generators):
print(f"generator {idx+1}: {generator}, paulix_op: {paulixops[idx]}")
paulix_sector = qml.qchem.optimal_sector(H, generators, n_electrons)
print(paulix_sector)
H_tapered = qml.taper(H, generators, paulixops, paulix_sector)
H_tapered_coeffs, H_tapered_ops = H_tapered.terms()
H_tapered = qml.Hamiltonian(jnp.real(jnp.array(H_tapered_coeffs)), H_tapered_ops)
print(H_tapered)
print (len(H_tapered))
H_sparse = qml.SparseHamiltonian(H.sparse_matrix(), wires=H.wires)
H_tapered_sparse = qml.SparseHamiltonian(H_tapered.sparse_matrix(), wires=H_tapered.wires)
print("Eigenvalues of H:\n", qml.eigvals(H_sparse, k=16))
print("\nEigenvalues of H_tapered:\n", qml.eigvals(H_tapered_sparse, k=4))
Tapering the reference state
state_tapered = qml.qchem.taper_hf(generators, paulixops, paulix_sector,
num_electrons=n_electrons, num_wires=len(H.wires))
print(state_tapered)
dev = qml.device("default.qubit", wires=H.wires)
@qml.qnode(dev, interface="jax")
def circuit():
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]), wires=H.wires)
return qml.state()
qubit_state = circuit()
HF_energy = qubit_state.T @ H.sparse_matrix().toarray() @ qubit_state
print(f"HF energy: {jnp.real(HF_energy):.8f} Ha")
dev = qml.device("lightning.qubit", wires=H_tapered.wires)
@qml.qnode(dev, interface="jax")
def circuit():
qml.BasisState(jnp.array([1, 1, 1, 1, 0, 0, 0, 0]), wires=H_tapered.wires)
return qml.state()
qubit_state = circuit()
HF_energy = qubit_state.T @ H_tapered.sparse_matrix().toarray() @ qubit_state
print(f"HF energy (tapered): {jnp.real(HF_energy):.8f} Ha")
VQE simulation
singles, doubles = qml.qchem.excitations(n_electrons, len(H.wires))
tapered_doubles = [
qml.taper_operation(qml.DoubleExcitation, generators, paulixops, paulix_sector,
wire_order=H.wires, op_wires=double) for double in doubles
]
tapered_singles = [
qml.taper_operation(qml.SingleExcitation, generators, paulixops, paulix_sector,
wire_order=H.wires, op_wires=single) for single in singles
]
dev = qml.device("lightning.qubit", wires=H_tapered.wires)
@qml.qnode(dev, interface="jax")
def tapered_circuit(params):
qml.BasisState(state_tapered, wires=H_tapered.wires)
for idx, tapered_op in enumerate(tapered_doubles + tapered_singles):
tapered_op(params[idx])
return qml.expval(H_tapered)
import optax
import catalyst
opt = optax.sgd(learning_rate=0.8) # sgd stands for StochasticGradientDescent
init_params = jnp.zeros(len(doubles) + len(singles))
def update_step(i, params, opt_state):
"""Perform a single gradient update step"""
grads = catalyst.grad(tapered_circuit)(params)
updates, opt_state = opt.update(grads, opt_state)
params = optax.apply_updates(params, updates)
return (params, opt_state)
loss_history = []
opt_state = opt.init(init_params)
params = init_params
for i in range(1, 101):
params, opt_state = update_step(i, params, opt_state)
energy = tapered_circuit(params)
if not i % 5:
print(f"n: {i}, E: {energy:.8f} Ha, Params: {params}")