This is the command to get the latest release of QuESTlink:

Import["https://qtechtheory.org/questlink.m"];

The next command downloads the pre-compiled single-thread simulator. For significantly faster multithreaded and GPU simulation, recompile QuESTlink (see guide here) and use CreateLocalQuESTEnv[]

CreateDownloadedQuESTEnv[];

Preparing an ansatz

Here’s how we could hard-code a 6 qubit, 12 parameter ansatz. The parameters are labelled θ[i] for convenience, but they could be any symbols we want (θi, (Subscript[θ, i]), etc)

Circuit[ 
	(Subscript[Rx, 0])[θ[1]](Subscript[Ry, 1])[θ[2]](Subscript[Ry, 2])[θ[3]](Subscript[Rx, 3])[θ[4]] (Subscript[Rx, 4])[θ[5]] (Subscript[Rx, 5])[θ[6]] 
	(Subscript[C, 0,1,2,4,5])[(Subscript[U, 3])[(({{Exp[I θ[7]], 0}, {0, Exp[-I θ[7]]}}))]]
	R[θ[7],(Subscript[Y, 0])(Subscript[X, 1])(Subscript[Z, 2])] R[θ[8], (Subscript[Y, 1])(Subscript[X, 2])(Subscript[Z, 3])]  R[θ[9], (Subscript[Y, 2])(Subscript[X, 3])(Subscript[Z, 4])] R[θ[10], (Subscript[Y, 3])(Subscript[X, 4])(Subscript[Z, 5])]
	(Subscript[C, 0,2,4])[R[θ[11], (Subscript[Z, 1])(Subscript[X, 3])(Subscript[Y, 5])]]
	R[θ[12], (Subscript[Y, 0])(Subscript[Y, 1])(Subscript[Y, 2])(Subscript[Y, 3])(Subscript[Y, 4])(Subscript[Y, 5])]];
DrawCircuit[%]

More often, we wish to generate an ansatz programmatically.

nQ = 6;
u = Block[{i=1}, Flatten @ {
Table[ (Subscript[g, q-1])[θ[i++]], {q,nQ}, {g,{Rx,Ry,Rz}}],
Table[(Subscript[C, DeleteCases[Range[0,nQ-1],q]])[(Subscript[U, q])[(({{Exp[I θ[i++]], 0}, {0, 1}})) ]], {q,0,nQ-1}],
Table[R[θ[i++], Product[(Subscript[σ, q-1]),{q,nQ}]],  {σ,{Y,X,Z}}]}];
nθ = Length[u];

DrawCircuit[u]

Since the parameters symbolic, we could study the ansatz operator analytically. Here’s its first element, when every second parameter has value π.

u  /. θ[_?EvenQ] -> π;
CalcCircuitMatrix[%]  //First//First

Preparing a Hamiltonian

We could hard-code our own Hamiltonian...

h = .1 (Subscript[X, 0])(Subscript[Y, 1])(Subscript[Z, 2]) + .3 (Subscript[Y, 0])(Subscript[X, 2]) - .4 (Subscript[Z, 0])(Subscript[Z, 1])(Subscript[Z, 2]);

CalcPauliSumMatrix[h] //Chop //MatrixForm

but more likely we’ll want to produce one from a file of coefficients and Pauli strings. Here we download a 369-term 6-qubit Lithium Hydride Hamiltonian

h = GetPauliSumFromCoeffs["https://questlink.qtechtheory.org/demo_hamiltonian.txt"];

Length[h]
h[[;;30]]

We can determine its ground-state exactly, which we’ll later compare to that which our variational algorithm estimates

gs = Min @ Eigenvalues @ CalcPauliSumMatrix[h]

Preparing quantum states

To perform numerical simulation of our ansatz , we need to create quantum registers, which are stored in the backend QuEST process.

ψ = CreateQureg[nQ];
InitPlusState[ψ];

GetQuregMatrix[ψ] // Chop

To apply our ansatz circuit to a state, we must replace the parameters with numerical values

ApplyCircuit[ u /. θ[_] :> RandomReal[], ψ];

GetQuregMatrix[ψ]

We must create a “working Qureg” hψ to compute expected values. We see our starting state isn’t especially close to the ground-state.

hψ = CreateQureg[nQ];
CalcExpecPauliSum[ψ, h, hψ]

We’ll also need additional Quregs to keep track of our fixed input state inψ = |+⟩, and nθ derivative states dψ[i] = (Ket[((∂ψ)/(∂   θ[i]))]) (where  ψ = u(θ) inψ )

inψ= CreateQureg[nQ];
InitPlusState[inψ];

dψ = CreateQuregs[nQ, nθ];

Computing variational observables

Let’s re-randomise the ansatz parameters.

vθ = Table[ θ[i]->RandomReal[], {i,nθ}]

We set each dψ[i] = (Ket[((∂ u)/(∂   θ[i])) inψ]) for the given assignment of θ

CalcQuregDerivs[u, inψ, vθ, dψ];

The imaginary time ‘tensor’ is simply Re[  ⟨((∂ψ)/(∂   θ[i]))||((∂ψ)/(∂   θ[j]))⟩ ]

m = Re @ CalcInnerProducts[dψ];
m //Chop //MatrixForm

We set |ψ⟩ = u(θ) |inψ⟩,  and |hψ⟩ = h |ψ⟩

ApplyCircuit[u /. vθ, CloneQureg[ψ,inψ]];
ApplyPauliSum[ψ, h, hψ];

The energy gradient is then simply Re[ ⟨ψ|h|((∂ψ)/(∂   θ[i]))⟩ ]

v = Re @ CalcInnerProducts[hψ, dψ];
v // MatrixForm

A single iteration of imaginary time evolution would use these observables to update the parameters via  m Δθ = -v Δt

Δθ = Δt LinearSolve[m, -v];
Δθ // MatrixForm

hence the updated parameters, dependent on Δt, would be

vθ[[All,2]] += Δθ

Performing variational imaginary time

We simply repeat these steps to perform variational imaginary time minimization. We’ll measure and record the energy at each iteration.

vθ = Table[ θ[i]->RandomReal[], {i,nθ}];

Δt = .01;
nt = 100;

en = Table[
	(* compute imaginary time tensor *)
	CalcQuregDerivs[u, inψ, vθ, dψ];
	m = Re @ CalcInnerProducts[dψ];

	(* compute energy gradient *)
	ApplyCircuit[u /. vθ, CloneQureg[ψ,inψ]];
	ApplyPauliSum[ψ, h, hψ];
	v = Re @ CalcInnerProducts[hψ, dψ];

	(* update parameters under imaginary time *)
	Δθ = Δt LinearSolve[m, -v];
	vθ[[All,2]] += Δθ;

	(* record the current energy *)
	CalcExpecPauliSum[ψ,h,hψ],

	(* perform nt iterations *)
	nt
];

Here’s how the energy evolved

ListLinePlot[en, 
	GridLines->{{},{{gs,Red}}}, 
	PlotRange->{gs,Automatic},
	AxesLabel-> {"Iteration", "Energy"}]

and the final energy, and its corresponding parameters

Row[{Last[en], "  /  ", gs}]
vθ