microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
samples/estimation/df-chemistry/src/DoubleFactorizedChemistry.qs
361lines · modecode
| 1 | // Copyright (c) Microsoft Corporation. All rights reserved. |
| 2 | // Licensed under the MIT License. |
| 3 | import Std.Arrays.*; |
| 4 | import Std.Canon.*; |
| 5 | import Std.Convert.*; |
| 6 | import Std.Diagnostics.*; |
| 7 | import Std.Math.*; |
| 8 | import Std.ResourceEstimation.*; |
| 9 | import Std.Arithmetic.*; |
| 10 | import Std.TableLookup.*; |
| 11 | import Prepare.*; |
| 12 | |
| 13 | // ------------------------------------------ // |
| 14 | // DF chemistry (public operations and types) // |
| 15 | // ------------------------------------------ // |
| 16 | |
| 17 | /// # Summary |
| 18 | /// The description of a double-factorized Hamiltonian |
| 19 | /// |
| 20 | /// # Reference |
| 21 | /// arXiv:2007.14460, p. 9, eq. 9 |
| 22 | struct DoubleFactorizedChemistryProblem { |
| 23 | // Number of orbitals (N, p. 8) |
| 24 | NumOrbitals : Int, |
| 25 | // one-body norm (ǁL⁽⁻¹⁾ǁ, p. 8, eq. 16) |
| 26 | OneBodyNorm : Double, |
| 27 | // one-body norm (¼∑ǁL⁽ʳ⁾ǁ², p. 8, eq. 16) |
| 28 | TwoBodyNorm : Double, |
| 29 | // eigenvalues in the EVD of the one-electron Hamiltonian (λₖ, p. 54, eq. 67) |
| 30 | OneBodyEigenValues : Double[], |
| 31 | // eigenvectors in the EVD of the one-electron Hamiltonian (Rₖ, p. 54, eq. 67) |
| 32 | OneBodyEigenVectors : Double[][], |
| 33 | // norms inside Λ_SH (p. 56, eq. 77) |
| 34 | Lambdas : Double[], |
| 35 | // eigenvalues in the EVDs of the two-electron Hamiltonian for all r (λₖ⁽ʳ⁾, p. 56, eq. 77) |
| 36 | TwoBodyEigenValues : Double[][], |
| 37 | // eigenvectors in the EVDs of the two-electron Hamiltonian for all r (R⁽ʳ⁾ₖ, p. 56, eq. 77) |
| 38 | TwoBodyEigenVectors : Double[][][], |
| 39 | } |
| 40 | |
| 41 | struct DoubleFactorizedChemistryParameters { |
| 42 | // Standard deviation (ΔE, p. 8, eq. 1) |
| 43 | // Typically set to 0.001 |
| 44 | StandardDeviation : Double, |
| 45 | } |
| 46 | |
| 47 | /// # Summary |
| 48 | /// Performs quantum circuit in (p. 45, eq. 33) for double-factorized |
| 49 | /// Hamiltonian; also prepares phase gradient registers to use phase |
| 50 | /// gradient technique (p. 55) |
| 51 | operation DoubleFactorizedChemistry( |
| 52 | problem : DoubleFactorizedChemistryProblem, |
| 53 | parameters : DoubleFactorizedChemistryParameters |
| 54 | ) : Unit { |
| 55 | let constants = ComputeConstants(problem, parameters); |
| 56 | |
| 57 | use register0 = Qubit[problem.NumOrbitals]; |
| 58 | use register1 = Qubit[problem.NumOrbitals]; |
| 59 | use phaseGradientRegister = Qubit[constants.RotationAngleBitPrecision]; |
| 60 | |
| 61 | // compute number of repetitions |
| 62 | let norm = problem.OneBodyNorm + problem.TwoBodyNorm; |
| 63 | let repetitions = Ceiling((PI() * norm) / (2.0 * parameters.StandardDeviation)); |
| 64 | |
| 65 | let walkStep = MakeWalkStep(problem, constants); |
| 66 | use walkStepHelper = Qubit[walkStep.NGarbageQubits]; |
| 67 | |
| 68 | within { |
| 69 | X(phaseGradientRegister[0]); |
| 70 | Adjoint ApplyQFT(Reversed(phaseGradientRegister)); |
| 71 | } apply { |
| 72 | within { |
| 73 | RepeatEstimates(repetitions); |
| 74 | } apply { |
| 75 | walkStep.StepOp(register0, register1, phaseGradientRegister, walkStepHelper); |
| 76 | } |
| 77 | } |
| 78 | } |
| 79 | |
| 80 | // -------------------------------------------- // |
| 81 | // DF chemistry (internal operations and types) // |
| 82 | // -------------------------------------------- // |
| 83 | |
| 84 | /// # Summary |
| 85 | /// Constants that are used in the implementation of the one- and two- |
| 86 | /// electron operators, which are computed from the double factorized |
| 87 | /// problem and parameters. |
| 88 | internal struct DoubleFactorizedChemistryConstants { |
| 89 | RotationAngleBitPrecision : Int, |
| 90 | StatePreparationBitPrecision : Int, |
| 91 | TargetError : Double |
| 92 | } |
| 93 | |
| 94 | internal function ComputeConstants( |
| 95 | problem : DoubleFactorizedChemistryProblem, |
| 96 | parameters : DoubleFactorizedChemistryParameters |
| 97 | ) : DoubleFactorizedChemistryConstants { |
| 98 | let pj = 0.1; |
| 99 | let barEpsilon = Sqrt(pj); |
| 100 | let norm = problem.OneBodyNorm + problem.TwoBodyNorm; |
| 101 | let epsilon = 0.1 * parameters.StandardDeviation / norm; |
| 102 | |
| 103 | let RotationAngleBitPrecision = Ceiling(1.152 + Lg(Sqrt((IntAsDouble((problem.NumOrbitals - 1) * 8) * PI() * norm) / parameters.StandardDeviation)) + 0.5 * Lg(1.0 / barEpsilon)); |
| 104 | let StatePreparationBitPrecision = Ceiling(Lg(1.0 / epsilon) + 2.5); |
| 105 | let TargetError = 2.0^IntAsDouble(1 - StatePreparationBitPrecision); |
| 106 | |
| 107 | new DoubleFactorizedChemistryConstants { |
| 108 | RotationAngleBitPrecision = RotationAngleBitPrecision, |
| 109 | StatePreparationBitPrecision = StatePreparationBitPrecision, |
| 110 | TargetError = TargetError |
| 111 | } |
| 112 | } |
| 113 | |
| 114 | internal struct WalkStep { |
| 115 | NGarbageQubits : Int, |
| 116 | StepOp : (Qubit[], Qubit[], Qubit[], Qubit[]) => Unit |
| 117 | } |
| 118 | |
| 119 | /// # Summary |
| 120 | /// Performs one B[H] and reflection (p. 45, eq. 33) |
| 121 | internal function MakeWalkStep( |
| 122 | problem : DoubleFactorizedChemistryProblem, |
| 123 | constants : DoubleFactorizedChemistryConstants |
| 124 | ) : WalkStep { |
| 125 | let oneElectronOperator = MakeOneElectronOperator(problem.OneBodyEigenValues, problem.OneBodyEigenVectors, constants); |
| 126 | let twoElectronOperator = MakeTwoElectronOperator(problem, constants); |
| 127 | |
| 128 | let NGarbageQubits = oneElectronOperator.NGarbageQubits + twoElectronOperator.NGarbageQubits; |
| 129 | |
| 130 | new WalkStep { |
| 131 | NGarbageQubits = NGarbageQubits, |
| 132 | StepOp = WalkStepOperation(problem, oneElectronOperator, twoElectronOperator, _, _, _, _) |
| 133 | } |
| 134 | } |
| 135 | |
| 136 | internal operation WalkStepOperation( |
| 137 | problem : DoubleFactorizedChemistryProblem, |
| 138 | oneElectronOperator : OneElectronOperator, |
| 139 | twoElectronOperator : TwoElectronOperator, |
| 140 | register0 : Qubit[], |
| 141 | register1 : Qubit[], |
| 142 | phaseGradientRegister : Qubit[], |
| 143 | helper : Qubit[] |
| 144 | ) : Unit { |
| 145 | use ctl = Qubit(); |
| 146 | |
| 147 | let helperParts = Partitioned([oneElectronOperator.NGarbageQubits, twoElectronOperator.NGarbageQubits], helper); |
| 148 | Fact(IsEmpty(Tail(helperParts)), "wrong number of helper qubits in WalkStepOperation"); |
| 149 | |
| 150 | // apply Hamiltionian |
| 151 | within { |
| 152 | PrepareSingleQubit(problem.OneBodyNorm, problem.TwoBodyNorm, ctl); |
| 153 | } apply { |
| 154 | within { X(ctl); } apply { |
| 155 | Controlled oneElectronOperator.Apply([ctl], (register0, register1, phaseGradientRegister, [], helperParts[0])); |
| 156 | } |
| 157 | Controlled twoElectronOperator.Apply([ctl], (register0, register1, phaseGradientRegister, helperParts[1])); |
| 158 | } |
| 159 | |
| 160 | // reflection |
| 161 | ReflectAboutInteger(0, helper); |
| 162 | } |
| 163 | |
| 164 | internal struct OneElectronOperator { |
| 165 | NGarbageQubits : Int, |
| 166 | Apply : (Qubit[], Qubit[], Qubit[], Qubit[], Qubit[]) => Unit is Adj + Ctl |
| 167 | } |
| 168 | |
| 169 | /// # Summary |
| 170 | /// Performs quantum circuit for one-electron operator (p. 54, eq. 70) |
| 171 | internal function MakeOneElectronOperator( |
| 172 | eigenValues : Double[], |
| 173 | eigenVectors : Double[][], |
| 174 | constants : DoubleFactorizedChemistryConstants |
| 175 | ) : OneElectronOperator { |
| 176 | let data = Mapped(eigenValue -> [eigenValue < 0.0], eigenValues); |
| 177 | let prepare = MakePrepareArbitrarySuperpositionWithData(constants.TargetError, eigenValues, data); |
| 178 | |
| 179 | new OneElectronOperator { |
| 180 | NGarbageQubits = prepare.NGarbageQubits, |
| 181 | Apply = OneElectronOperatorOperation(eigenVectors, constants, prepare, _, _, _, _, _) |
| 182 | } |
| 183 | } |
| 184 | |
| 185 | internal operation OneElectronOperatorOperation( |
| 186 | eigenVectors : Double[][], |
| 187 | constants : DoubleFactorizedChemistryConstants, |
| 188 | prepare : PrepareArbitrarySuperposition, |
| 189 | register0 : Qubit[], |
| 190 | register1 : Qubit[], |
| 191 | phaseGradientRegister : Qubit[], |
| 192 | offsets : Qubit[], |
| 193 | helper : Qubit[] |
| 194 | ) : Unit is Adj + Ctl { |
| 195 | if BeginEstimateCaching("OneElectronOperator", IsEmpty(offsets) ? 0 | 1) { |
| 196 | // assertions |
| 197 | Fact(Length(helper) == prepare.NGarbageQubits, "invalid number of helper qubits in OneElectronOperatorOperation"); |
| 198 | |
| 199 | let precision = constants.RotationAngleBitPrecision; |
| 200 | let bitstrings = AllEigenVectorsAsBitString(eigenVectors, precision); |
| 201 | |
| 202 | use prepQubits = Qubit[prepare.NIndexQubits]; |
| 203 | use rotationQubits = Qubit[Length(Head(bitstrings))]; |
| 204 | use sign = Qubit(); |
| 205 | use spin = Qubit(); |
| 206 | |
| 207 | within { |
| 208 | if not IsEmpty(offsets) { |
| 209 | prepare.PrepareWithSelect(SelectWithOffset(_, offsets, _, _), prepQubits, [sign], helper); |
| 210 | } else { |
| 211 | prepare.Prepare(prepQubits, [sign], helper); |
| 212 | } |
| 213 | H(spin); |
| 214 | for i in IndexRange(register0) { |
| 215 | Controlled SWAP([spin], (register0[i], register1[i])); |
| 216 | } |
| 217 | if not IsEmpty(offsets) { |
| 218 | RippleCarryCGIncByLE(offsets, prepQubits); |
| 219 | } |
| 220 | Select(bitstrings, prepQubits, rotationQubits); |
| 221 | } apply { |
| 222 | ApplyGivensRotations(phaseGradientRegister, Chunks(precision, rotationQubits), register0); |
| 223 | Z(sign); |
| 224 | } |
| 225 | |
| 226 | EndEstimateCaching(); |
| 227 | } |
| 228 | } |
| 229 | |
| 230 | internal operation SelectWithOffset(data : Bool[][], offset : Qubit[], address : Qubit[], target : Qubit[]) : Unit is Adj + Ctl { |
| 231 | within { |
| 232 | RippleCarryCGIncByLE(offset, address); |
| 233 | } apply { |
| 234 | Select(data, address, target); |
| 235 | } |
| 236 | } |
| 237 | |
| 238 | internal struct TwoElectronOperator { |
| 239 | NGarbageQubits : Int, |
| 240 | Apply : (Qubit[], Qubit[], Qubit[], Qubit[]) => Unit is Ctl |
| 241 | } |
| 242 | |
| 243 | /// # Summary |
| 244 | /// Performs quantum circuit for two-electron operator (p. 57, eq. 78) |
| 245 | internal function MakeTwoElectronOperator(problem : DoubleFactorizedChemistryProblem, constants : DoubleFactorizedChemistryConstants) : TwoElectronOperator { |
| 246 | let lambdaPrepare = MakePrepareArbitrarySuperposition( |
| 247 | constants.TargetError, |
| 248 | problem.Lambdas |
| 249 | ); |
| 250 | |
| 251 | let eigenValuesFlattened = Flattened(problem.TwoBodyEigenValues); |
| 252 | let eigenVectorsFlattened = Flattened(problem.TwoBodyEigenVectors); |
| 253 | let oneElectronOperator = MakeOneElectronOperator(eigenValuesFlattened, eigenVectorsFlattened, constants); |
| 254 | |
| 255 | let numGarbageQubits = lambdaPrepare.NGarbageQubits + oneElectronOperator.NGarbageQubits; |
| 256 | |
| 257 | new TwoElectronOperator { NGarbageQubits = numGarbageQubits, Apply = TwoElectronOperatorOperation(problem, lambdaPrepare, oneElectronOperator, _, _, _, _) } |
| 258 | } |
| 259 | |
| 260 | internal operation TwoElectronOperatorOperation( |
| 261 | problem : DoubleFactorizedChemistryProblem, |
| 262 | lambdaPrepare : PrepareArbitrarySuperposition, |
| 263 | oneElectronOperator : OneElectronOperator, |
| 264 | register0 : Qubit[], |
| 265 | register1 : Qubit[], |
| 266 | phaseGradientRegister : Qubit[], |
| 267 | helper : Qubit[] |
| 268 | ) : Unit is Ctl { |
| 269 | let helperParts = Partitioned([lambdaPrepare.NGarbageQubits], helper); |
| 270 | let dataOffsets = ComputeOffsetDataSet(Mapped(Length, problem.TwoBodyEigenValues)); |
| 271 | |
| 272 | use rankQubits = Qubit[lambdaPrepare.NIndexQubits]; |
| 273 | use offsetQubits = Qubit[Length(dataOffsets[0])]; |
| 274 | |
| 275 | within { |
| 276 | lambdaPrepare.Prepare(rankQubits, [], helperParts[0]); |
| 277 | Select(dataOffsets, rankQubits, offsetQubits); |
| 278 | oneElectronOperator.Apply(register0, register1, phaseGradientRegister, offsetQubits, helperParts[1]); |
| 279 | } apply { |
| 280 | ReflectAboutInteger(0, helperParts[1]); |
| 281 | } |
| 282 | } |
| 283 | |
| 284 | internal function ComputeOffsetDataSet(lengths : Int[]) : Bool[][] { |
| 285 | mutable currentOffset = 0; |
| 286 | mutable offsets = [0, size = Length(lengths)]; |
| 287 | |
| 288 | for i in IndexRange(lengths) { |
| 289 | set offsets w/= i <- currentOffset; |
| 290 | set currentOffset += lengths[i]; |
| 291 | } |
| 292 | |
| 293 | let offsetWidth = Ceiling(Lg(IntAsDouble(currentOffset))); |
| 294 | return Mapped(IntAsBoolArray(_, offsetWidth), offsets); |
| 295 | } |
| 296 | |
| 297 | /// # Summary |
| 298 | /// Performs quantum circuit for phase rotations (p. 54, eq. 73) |
| 299 | internal operation ApplyGivensRotations(pgr : Qubit[], rotationQubits : Qubit[][], target : Qubit[]) : Unit is Adj + Ctl { |
| 300 | within { |
| 301 | let windows = Windows(2, target); |
| 302 | for i in RangeReverse(IndexRange(rotationQubits)) { |
| 303 | let rotationControls = rotationQubits[i]; |
| 304 | let qs = windows[i]; |
| 305 | ApplyWithMajoranaCliffords(1, ApplyRotationUsingRippleCarryAddition(pgr, rotationControls, _), qs); |
| 306 | ApplyWithMajoranaCliffords(0, Adjoint ApplyRotationUsingRippleCarryAddition(pgr, rotationControls, _), qs); |
| 307 | } |
| 308 | } apply { |
| 309 | X(Head(target)); |
| 310 | Y(Head(target)); |
| 311 | } |
| 312 | } |
| 313 | |
| 314 | internal operation ApplyRotationUsingRippleCarryAddition(pgr : Qubit[], rotationQubits : Qubit[], target : Qubit) : Unit is Adj { |
| 315 | Controlled Adjoint RippleCarryCGIncByLE([target], (Reversed(Rest(rotationQubits)), pgr)); |
| 316 | } |
| 317 | |
| 318 | internal operation ApplyWithMajoranaCliffords(x : Int, op : (Qubit => Unit is Adj), qs : Qubit[]) : Unit is Adj { |
| 319 | // x is the Majorana type (which can be either 0 or 1, see p. 52, eq. 59) |
| 320 | within { |
| 321 | Adjoint S(qs[x]); |
| 322 | H(qs[x]); |
| 323 | H(qs[1 - x]); |
| 324 | CNOT(qs[1], qs[0]); |
| 325 | } apply { |
| 326 | op(qs[0]); |
| 327 | } |
| 328 | } |
| 329 | |
| 330 | internal function AllEigenVectorsAsBitString(eigenVectors : Double[][], precision : Int) : Bool[][] { |
| 331 | mutable bitstrings = []; |
| 332 | |
| 333 | let tau = 2.0 * PI(); |
| 334 | let preFactor = 2.0^IntAsDouble(precision); |
| 335 | |
| 336 | for eigenVector in eigenVectors { |
| 337 | // Computes rotation angles for Majorana operator ($\vec u$ in p. 52, eq. 55) |
| 338 | mutable result = []; |
| 339 | mutable sins = 1.0; |
| 340 | |
| 341 | for index in 0..Length(eigenVector) - 2 { |
| 342 | // We apply MinD, such that rounding errors do not lead to |
| 343 | // an argument for ArcCos which is larger than 1.0. (p. 52, eq. 56) |
| 344 | let theta = not (AbsD(sins) > 0.0) ? 0.0 | 0.5 * ArcCos(MinD(eigenVector[index] / sins, 1.0)); |
| 345 | |
| 346 | // all angles as bit string |
| 347 | let factor = theta / tau; |
| 348 | set result += Reversed(IntAsBoolArray(IsNaN(factor) ? 0 | Floor(preFactor * factor), precision)); |
| 349 | |
| 350 | set sins *= Sin(2.0 * theta); |
| 351 | } |
| 352 | |
| 353 | set bitstrings += [result]; |
| 354 | } |
| 355 | |
| 356 | bitstrings |
| 357 | } |
| 358 | |
| 359 | internal function IsNaN(value : Double) : Bool { |
| 360 | not (value == value) |
| 361 | } |
| 362 | |