# Copyright 2021 Huawei Technologies Co., Ltd
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ============================================================================
"""UCCSD ansatz."""
from collections import OrderedDict as ordict
import itertools
import numpy as np
from openfermion.chem import MolecularData
from openfermion.ops import FermionOperator
from openfermion.utils.indexing import down_index
from openfermion.utils.indexing import up_index
from openfermion.transforms import jordan_wigner
from mindquantum.gate import RX
from mindquantum.gate import H
from mindquantum.gate import X
from mindquantum.gate import RZ
from mindquantum.gate import RY
from .circuit import Circuit
def _para_uccsd_singlet_generator(mol, th=0):
"""
Generate a uccsd quantum circuit.
Args:
mol (molecular): A hdf5 molecular file generated by HiQ Fermion.
th (int, optional): A threadhold of parameters. If a parameter is
lower than the threadhold, than we will not update it in VQE
algorithm. Default: 0.
"""
n_qubits = mol.n_qubits
n_electrons = mol.n_electrons
params = {}
if n_qubits % 2 != 0:
raise ValueError('The total number of spin-orbitals should be even.')
out = []
out_tmp = []
n_spatial_orbitals = n_qubits // 2
n_occupied = int(np.ceil(n_electrons / 2))
n_virtual = n_spatial_orbitals - n_occupied
# Unpack amplitudes
n_single_amplitudes = n_occupied * n_virtual
# Generate excitations
spin_index_functions = [up_index, down_index]
# Generate all spin-conserving single and double excitations derived
# from one spatial occupied-virtual pair
for i, (p, q) in enumerate(
itertools.product(range(n_virtual), range(n_occupied))):
# Get indices of spatial orbitals
virtual_spatial = n_occupied + p
occupied_spatial = q
virtual_up = virtual_spatial * 2
occupied_up = occupied_spatial * 2
virtual_down = virtual_spatial * 2 + 1
occupied_down = occupied_spatial * 2 + 1
single_amps = mol.ccsd_single_amps[virtual_up, occupied_up]
double1_amps = mol.ccsd_double_amps[virtual_up, occupied_up,
virtual_down, occupied_down]
single_amps_name = 'p' + str(i)
double1_amps_name = 'p' + str(i + n_single_amplitudes)
for spin in range(2):
# Get the functions which map a spatial orbital index to a
# spin orbital index
this_index = spin_index_functions[spin]
other_index = spin_index_functions[1 - spin]
# Get indices of spin orbitals
virtual_this = this_index(virtual_spatial)
virtual_other = other_index(virtual_spatial)
occupied_this = this_index(occupied_spatial)
occupied_other = other_index(occupied_spatial)
# Generate single excitations
if abs(single_amps) > th:
params[single_amps_name] = single_amps
fermion_ops1 = FermionOperator(
((occupied_this, 1), (virtual_this, 0)), 1)
fermion_ops2 = FermionOperator(
((virtual_this, 1), (occupied_this, 0)), 1)
out.append([fermion_ops1 - fermion_ops2, single_amps_name])
# Generate double excitation
if abs(double1_amps) > th:
params[double1_amps_name] = double1_amps
fermion_ops1 = FermionOperator(
((virtual_this, 1), (occupied_this, 0), (virtual_other, 1),
(occupied_other, 0)), 1)
fermion_ops2 = FermionOperator(
((occupied_other, 1), (virtual_other, 0),
(occupied_this, 1), (virtual_this, 0)), 1)
out.append([fermion_ops1 - fermion_ops2, double1_amps_name])
out.extend(out_tmp)
out_tmp = []
# Generate all spin-conserving double excitations derived
# from two spatial occupied-virtual pairs
for i, ((p, q), (r, s)) in enumerate(
itertools.combinations(
itertools.product(range(n_virtual), range(n_occupied)), 2)):
# Get indices of spatial orbitals
virtual_spatial_1 = n_occupied + p
occupied_spatial_1 = q
virtual_spatial_2 = n_occupied + r
occupied_spatial_2 = s
virtual_1_up = virtual_spatial_1 * 2
occupied_1_up = occupied_spatial_1 * 2
virtual_2_up = virtual_spatial_2 * 2 + 1
occupied_2_up = occupied_spatial_2 * 2 + 1
double2_amps = mol.ccsd_double_amps[virtual_1_up, occupied_1_up,
virtual_2_up, occupied_2_up]
double2_amps_name = 'p' + str(i + 2 * n_single_amplitudes)
# Generate double excitations
for (spin_a, spin_b) in itertools.product(range(2), repeat=2):
# Get the functions which map a spatial orbital index to a
# spin orbital index
index_a = spin_index_functions[spin_a]
index_b = spin_index_functions[spin_b]
# Get indices of spin orbitals
virtual_1_a = index_a(virtual_spatial_1)
occupied_1_a = index_a(occupied_spatial_1)
virtual_2_b = index_b(virtual_spatial_2)
occupied_2_b = index_b(occupied_spatial_2)
if abs(double2_amps) > th:
params[double2_amps_name] = double2_amps
fermion_ops1 = FermionOperator(
((virtual_1_a, 1), (occupied_1_a, 0), (virtual_2_b, 1),
(occupied_2_b, 0)), 1)
fermion_ops2 = FermionOperator(
((occupied_2_b, 1), (virtual_2_b, 0), (occupied_1_a, 1),
(virtual_1_a, 0)), 1)
out.append([fermion_ops1 - fermion_ops2, double2_amps_name])
return out, params
def _transform2pauli(fermion_ansatz):
"""
Transform a fermion ansatz to pauli ansatz based on jordan-wigner
transformation.
"""
out = ordict()
for i in fermion_ansatz:
qubit_generator = jordan_wigner(i[0])
if qubit_generator.terms != {}:
for key, term in qubit_generator.terms.items():
if key not in out:
out[key] = ordict({i[1]: float(term.imag)})
else:
if i[1] in out[key]:
out[key][i[1]] += float(term.imag)
else:
out[key][i[1]] = float(term.imag)
return out
def _pauli2circuit(pauli_ansatz):
"""Transform a pauli ansatz to parameterized quantum circuit."""
circuit = Circuit()
for k, v in pauli_ansatz.items():
circuit += decompose_single_term_time_evolution(k, v)
return circuit
[docs]def decompose_single_term_time_evolution(term, para):
"""
Decompose a time evolution gate into basic quantum gates.
This function only work for the hamiltonian with only single pauli word.
For example, exp(-i * t * ham), ham can only be a single pauli word, such
as ham = X0 x Y1 x Z2, and at this time, term will be
((0, 'X'), (1, 'Y'), (2, 'Z')). When the evolution time is expressd as
t = a*x + b*y, para would be {'x':a, 'y':b}.
Args:
term (tuple, QubitOperator): the hamiltonian term of just the
evolution qubit operator.
para (dict): the parameters of evolution operator.
Returns:
Circuit, a quantum circuit.
Example:
>>> from projectq.ops import QubitOperator
>>> ham = QubitOperator('X0 Y1')
>>> circuit = decompose_single_term_time_evolution(ham, {'a':1})
>>> print(circuit)
H(0)
RX(1.571,1)
X(1 <-: 0)
RZ(a|1)
X(1 <-: 0)
RX(10.996,1)
H(0)
"""
if not isinstance(term, tuple):
try:
if len(term.terms) != 1:
raise ValueError("Only work for single term time \
evolution operator, but get {}".format(len(term)))
term = list(term.terms.keys())[0]
except TypeError:
raise Exception("Not supported type:{}".format(type(term)))
out = []
term = sorted(term)
rxs = []
if len(term) == 1: # single pauli operator
if term[0][1] == 'X':
out.append(RX(para).on(term[0][0]))
elif term[0][1] == 'Y':
out.append(RY(para).on(term[0][0]))
else:
out.append(RZ(para).on(term[0][0]))
else:
for index, action in term:
if action == 'X':
out.append(H.on(index))
elif action == 'Y':
rxs.append(len(out))
out.append(RX(np.pi / 2).on(index))
for i in range(len(term) - 1):
out.append(X.on(term[i + 1][0], term[i][0]))
out.append(RZ({i: 2 * j for i, j in para.items()}).on(term[-1][0]))
for i in range(len(out) - 1)[::-1]:
if i in rxs:
out.append(RX(np.pi * 3.5).on(out[i].obj_qubits))
else:
out.append(out[i])
return Circuit(out)
[docs]def generate_uccsd(filename, th=0):
"""
Generate a uccsd quantum circuit based on a molecular data generated by
HiQfermion or openfermion.
Args:
filename (str): the name of the molecular data file.
th (int): the threshold to filt the uccsd amplitude. When th < 0, we
will keep all amplitudes. When th == 0, we will keep all amplitude
that are positive.
Returns:
- **uccsd_circuit** (Circuit), the ansatz circuit generated by uccsd method.
- **initial_amplitudes** (numpy.ndarray), the initial parameter values of uccsd circuit.
- **parameters_name** (list[str]), the name of initial parameters.
- **qubit_hamiltonian** (QubitOperator), the hamiltonian of the molecule.
- **n_qubits** (int), the number of qubits in simulation.
- **n_electrons**, the number of electrons of the molecule.
"""
mol = MolecularData(filename=filename)
mol.load()
print("ccsd:{}.".format(mol.ccsd_energy))
print("fci:{}.".format(mol.fci_energy))
fermion_ansatz, parameters = _para_uccsd_singlet_generator(mol, th)
pauli_ansatz = _transform2pauli(fermion_ansatz)
uccsd_circuit = _pauli2circuit(pauli_ansatz)
qubit_hamiltonian = jordan_wigner(mol.get_molecular_hamiltonian())
qubit_hamiltonian.compress()
parameters_name = list(parameters.keys())
initial_amplitudes = [parameters[i] for i in parameters_name]
return uccsd_circuit, \
initial_amplitudes, \
parameters_name, \
qubit_hamiltonian, \
mol.n_qubits, \
mol.n_electrons