Skip to content

Commit

Permalink
Merge pull request #27 from yipenghuang0302/master
Browse files Browse the repository at this point in the history 
Checking in code for validated benchmark programs: H2 ground state estimation, Shor's, Grover's
  • Loading branch information
EPiQC authored Jul 8, 2019
2 parents 6ff325d + 0a98828 commit 98275a9
Show file tree
Hide file tree
Showing 26 changed files with 1,643 additions and 299 deletions.
157 changes: 76 additions & 81 deletions Algorithms/Ground_State_Estimation/ground_state_estimation.h2.scaffold
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,12 @@
#include <math.h>
#include <../../Library/cRz/cRz.scaffold>

#define M 4 //number of bottom register qubits
#define TROTTER 16 //number of bottom register qubits
#define TROTTER 8

#define _upper_sz 4
#define _iteration 0
#define _estimate 0

enum spin {alpha, beta};
const enum spin sigma[M] = {
Expand All @@ -16,54 +21,57 @@ const enum spin sigma[M] = {
// info=0.00,0.7414
// nuc=0.713776188
// Ehf=-1.116685636
// 0,0=-1.252477495
// 1,1=-0.475934275
// 0,0,0,0=0.674493166
// 0,1,0,1=0.181287518
// 0,1,1,0=0.663472101
// 0,0=-1.252477495 -1.2563390710389655
// 1,1=-0.475934275 -0.4718960093502949
// 0,0,0,0=0.674493166 0.6757101541858549
// 0,1,0,1=0.181287518 0.18093119996471013
// 0,1,1,0=0.663472101 0.664581729691277
// 1,1,1,1=0.697398010

// h_double[0][1][1][0]=h_double[1][0][0][1]=0.674493166;
// h_double[0][2][2][0]=h_double[2][0][0][2]=0.663472101;
// h_double[0][3][3][0]=h_double[3][0][0][3]=0.663472101;
// h_double[1][2][2][1]=h_double[2][1][1][2]=0.663472101;
// h_double[1][3][3][1]=h_double[3][1][1][3]=0.663472101;
// h_double[0][1][1][0]=h_double[1][0][0][1]=0.6757101541858549;
// h_double[0][2][2][0]=h_double[2][0][0][2]=0.664581729691277;
// h_double[0][3][3][0]=h_double[3][0][0][3]=0.664581729691277;
// h_double[1][2][2][1]=h_double[2][1][1][2]=0.664581729691277;
// h_double[1][3][3][1]=h_double[3][1][1][3]=0.664581729691277;
// h_double[2][3][3][2]=h_double[3][2][2][3]=0.697398010;
//
// h_double[0][2][0][2]=h_double[1][3][1][3]=0.181287518;
// h_double[2][0][2][0]=h_double[3][1][3][1]=0.000000;
// h_double[0][3][1][2]=h_double[0][1][3][2]=0.181287518;
// h_double[0][2][0][2]=h_double[1][3][1][3]=0.18093119996471013;
// h_double[2][0][2][0]=h_double[3][1][3][1]=0.18093119996471013;
// h_double[0][3][1][2]=h_double[0][1][3][2]=0.18093119996471013;

// according to 1.4768229
// h_double[2][1][3][0]=h_double[2][3][1][0]=0.18093119996471013;

const double h_single[M][M] = {
{-1.252477495, 0.0, 0.0, 0.0},
{0.0, -1.252477495, 0.0, 0.0},
{0.0, 0.0, -0.475934275, 0.0},
{0.0, 0.0, 0.0, -0.475934275}
{-1.2563390710389655, 0.0, 0.0, 0.0},
{0.0, -1.2563390710389655, 0.0, 0.0},
{0.0, 0.0, -0.4718960093502949, 0.0},
{0.0, 0.0, 0.0, -0.4718960093502949}
};

const double h_double[M][M][M][M] =
{
{
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.674493166,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.181287518,0.0}},
{{0.0,0.0,0.181287518,0.0},{0.0,0.0,0.0,0.0},{0.663472101,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.181287518,0.0},{0.0,0.0,0.0,0.0},{0.663472101,0.0,0.0,0.0}}
{{0.0,0.0,0.0,0.0},{0.6757101541858549,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.18093119996471013,0.0}},
{{0.0,0.0,0.18093119996471013,0.0},{0.0,0.0,0.0,0.0},{0.664581729691277,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.18093119996471013,0.0},{0.0,0.0,0.0,0.0},{0.664581729691277,0.0,0.0,0.0}}
},
{
{{0.0,0.674493166,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.6757101541858549,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.663472101,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.181287518},{0.0,0.0,0.0,0.0},{0.0,0.663472101,0.0,0.0}}
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.664581729691277,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.18093119996471013},{0.0,0.0,0.0,0.0},{0.0,0.664581729691277,0.0,0.0}}
},
{
{{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.663472101,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.664581729691277,0.0},{0.0,0.0,0.0,0.0},{0.18093119996471013,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.664581729691277,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.697398010,0.0}}
},
{
{{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.663472101},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.664581729691277},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.664581729691277},{0.0,0.0,0.0,0.0},{0.0,0.18093119996471013,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.697398010},{0.0,0.0,0.0,0.0}},
{{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0},{0.0,0.0,0.0,0.0}}
}
Expand All @@ -77,21 +85,7 @@ module ControlledPhase (
qbit target,
double theta
) {
Rz(control, theta);
}

// Implements a controlled rotation
// inputs are the target qubit, control qubit and the rotation angle
module ControlledRotation (
qbit control,
qbit target,
double theta
) {
// Rz (control, theta / 2.0);
CNOT (control, target);
Rz (target, theta / 2.0);
CNOT (control, target);
Rz (target, -1.0 * theta / 2.0);
// Rz(control, theta);
}

// Implement the Y gate defined in the GFIs Y=Rx(-PI/2), we'll name it Y_Rx
Expand Down Expand Up @@ -139,10 +133,10 @@ module DoubleExcitationOperator (
}
CNOT(reg[p], reg[q]); CNOT(reg[q], reg[r]); CNOT(reg[r], reg[s]);
switch (i) { // M operator (as defined in A1)
case 0: ControlledRotation(control, reg[s], theta0); break;
case 1: ControlledRotation(control, reg[s], theta1); break;
case 2: ControlledRotation(control, reg[s], theta2); break;
case 3: ControlledRotation(control, reg[s], theta3); break;
case 0: cRz(control, reg[s], theta0); break;
case 1: cRz(control, reg[s], theta1); break;
case 2: cRz(control, reg[s], theta2); break;
case 3: cRz(control, reg[s], theta3); break;
}
CNOT(reg[r], reg[s]); CNOT(reg[q], reg[r]); CNOT(reg[p], reg[q]);
switch (i) { // M^dagger operator is identical since matrices are Hermitian
Expand All @@ -158,19 +152,20 @@ int main ()
{

qbit bottomregister[M];
// for ( int i=0; i<M; i++ ) {
// PrepZ ( bottomregister[i], 1 );
// }
PrepZ ( bottomregister[3], 1 );
PrepZ ( bottomregister[2], 1 );

// eigenstate electron orbital assignments
// bonding orbitals
PrepZ ( bottomregister[0], 1 );
PrepZ ( bottomregister[1], 1 );
PrepZ ( bottomregister[0], 0 );
// anitbonding orbitals
PrepZ ( bottomregister[2], 0 );
PrepZ ( bottomregister[3], 0 );

qbit topregister[2];
PrepZ ( topregister[0], 0 );
H(topregister[0]);

for (int unitary=0; unitary<1; unitary++) {
for (int unitary=0; unitary<(1<<(_upper_sz-_iteration-1)); unitary++) {
double time_step = 1.0/(double)TROTTER;
for (int rep=0; rep<TROTTER; rep++) {

Expand All @@ -179,10 +174,9 @@ int main ()
ControlledPhase (
topregister[0],
bottomregister[i],
h_single[i][i] * time_step
0.5 * h_single[i][i] * time_step
);
}

// two electron operators: number-number operators
double theta = 0.0;
for (int q=0; q<M; q++) {
Expand All @@ -193,37 +187,43 @@ int main ()
}
}
}
Rz (topregister[0], theta * time_step); // FIXME: should really be phase

Rz (topregister[0], theta * time_step);

// single electron operators
for (int i=0; i<M; i++) {
cRz (
topregister[0],
bottomregister[i],
- 0.5 * h_single[i][i] * time_step
);
}

for (int p=0; p<M; p++) {
double thetas = 0.0;
for (int q=0; q<M; q++) {
if (p!=q) { // TODO: check this against UCSB documentation
// if (p<q) { // TODO: check this against UCSB documentation
thetas += h_double[p][q][q][p]; // TODO: this might be accidentally too many
if (p!=q) { // DONE: check this against UCSB documentation
thetas += h_double[p][q][q][p];
if (sigma[p]==sigma[q]) {
thetas -= h_double[p][q][p][q]; // TODO: h2020 and h3131 are missing
thetas -= h_double[p][q][p][q]; // h2020 and h3131 are missing; McArdle confirms they are zero
}
}
}

ControlledRotation (
cRz (
topregister[0],
bottomregister[p],
-1.0 * thetas * time_step
- 0.25 * thetas * time_step
);
}


// Cross validated against Seeley 2012
for (int t=3; t>=0; t--) {
for (int c=t-1; c>=0; c--) {
CNOT(bottomregister[c], bottomregister[t]);

double eta = 0.25 * h_double[c][t][t][c];
if (sigma[c]==sigma[t]) eta -= 0.25*h_double[c][t][c][t];
ControlledRotation(topregister[0], bottomregister[t], 2*eta*time_step);
cRz(topregister[0], bottomregister[t], eta*time_step);

CNOT(bottomregister[c], bottomregister[t]);
}
Expand All @@ -234,24 +234,19 @@ int main ()
topregister[0],
bottomregister,
0, 1, 2, 3,
-time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0,
-time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0,
time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0,
time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 4.0
-time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0,
-time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0,
time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0,
time_step * (h_double[0][3][1][2]+h_double[0][1][3][2]) / 8.0
);
}
}

// Rz(topregister[0],-M_PI/64);//
// Rz(topregister[0],-M_PI/32);//
// Rz(topregister[0],-M_PI/16);//
// Rz(topregister[0],-M_PI/8);//
// Rz(topregister[0],-M_PI/4);//
// Rz(topregister[0],-M_PI/2);//
for (int c=0; c<_iteration; c++) {
if ((_estimate>>c)&1) {
Rz ( topregister[0], -M_PI/pow(2,_iteration-c) );
}
}
H(topregister[0]);
// MeasZ(bottomregister[0]);
// MeasZ(bottomregister[1]);
// MeasZ(bottomregister[2]);
// MeasZ(bottomregister[3]);
// MeasZ(topregister[0]);

}
Loading

0 comments on commit 98275a9

Please sign in to comment.