microsoft/qdk
Publicmirrored from https://github.com/microsoft/qdkAvailable
samples/estimation/df-chemistry/src/df_chemistry.qs
361lines · modecode
| 1 | // Copyright (c) Microsoft Corporation. All rights reserved. |
| 2 | // Licensed under the MIT License. |
| 3 | namespace Microsoft.Quantum.Applications.Chemistry { |
| 4 | open Microsoft.Quantum.Arrays; |
| 5 | open Microsoft.Quantum.Canon; |
| 6 | open Microsoft.Quantum.Convert; |
| 7 | open Microsoft.Quantum.Diagnostics; |
| 8 | open Microsoft.Quantum.Math; |
| 9 | open Microsoft.Quantum.ResourceEstimation; |
| 10 | open Microsoft.Quantum.Unstable.Arithmetic; |
| 11 | open Microsoft.Quantum.Unstable.TableLookup; |
| 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 | newtype 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 | newtype 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 newtype 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 |
| 104 | + Lg(Sqrt((IntAsDouble((problem::NumOrbitals - 1) * 8) * PI() * norm) / parameters::StandardDeviation)) |
| 105 | + 0.5 * Lg(1.0 / barEpsilon)); |
| 106 | let StatePreparationBitPrecision = Ceiling(Lg(1.0 / epsilon) + 2.5); |
| 107 | let TargetError = 2.0^IntAsDouble(1 - StatePreparationBitPrecision); |
| 108 | |
| 109 | DoubleFactorizedChemistryConstants( |
| 110 | RotationAngleBitPrecision, |
| 111 | StatePreparationBitPrecision, |
| 112 | TargetError |
| 113 | ) |
| 114 | } |
| 115 | |
| 116 | internal newtype WalkStep = ( |
| 117 | NGarbageQubits: Int, |
| 118 | StepOp: (Qubit[], Qubit[], Qubit[], Qubit[]) => Unit |
| 119 | ); |
| 120 | |
| 121 | /// # Summary |
| 122 | /// Performs one B[H] and reflection (p. 45, eq. 33) |
| 123 | internal function MakeWalkStep( |
| 124 | problem : DoubleFactorizedChemistryProblem, |
| 125 | constants : DoubleFactorizedChemistryConstants |
| 126 | ) : WalkStep { |
| 127 | let oneElectronOperator = MakeOneElectronOperator(problem::OneBodyEigenValues, problem::OneBodyEigenVectors, constants); |
| 128 | let twoElectronOperator = MakeTwoElectronOperator(problem, constants); |
| 129 | |
| 130 | let NGarbageQubits = oneElectronOperator::NGarbageQubits + twoElectronOperator::NGarbageQubits; |
| 131 | |
| 132 | WalkStep(NGarbageQubits, WalkStepOperation(problem, oneElectronOperator, twoElectronOperator, _, _, _, _)) |
| 133 | } |
| 134 | |
| 135 | internal operation WalkStepOperation( |
| 136 | problem : DoubleFactorizedChemistryProblem, |
| 137 | oneElectronOperator : OneElectronOperator, |
| 138 | twoElectronOperator : TwoElectronOperator, |
| 139 | register0 : Qubit[], |
| 140 | register1 : Qubit[], |
| 141 | phaseGradientRegister : Qubit[], |
| 142 | helper : Qubit[] |
| 143 | ) : Unit { |
| 144 | use ctl = Qubit(); |
| 145 | |
| 146 | let helperParts = Partitioned([oneElectronOperator::NGarbageQubits, twoElectronOperator::NGarbageQubits], helper); |
| 147 | Fact(IsEmpty(Tail(helperParts)), "wrong number of helper qubits in WalkStepOperation"); |
| 148 | |
| 149 | // apply Hamiltionian |
| 150 | within { |
| 151 | PrepareSingleQubit(problem::OneBodyNorm, problem::TwoBodyNorm, ctl); |
| 152 | } apply { |
| 153 | within { X(ctl); } apply { |
| 154 | Controlled oneElectronOperator::Apply([ctl], (register0, register1, phaseGradientRegister, [], helperParts[0])); |
| 155 | } |
| 156 | Controlled twoElectronOperator::Apply([ctl], (register0, register1, phaseGradientRegister, helperParts[1])); |
| 157 | } |
| 158 | |
| 159 | // reflection |
| 160 | ReflectAboutInteger(0, helper); |
| 161 | } |
| 162 | |
| 163 | internal newtype OneElectronOperator = ( |
| 164 | NGarbageQubits: Int, |
| 165 | Apply: (Qubit[], Qubit[], Qubit[], Qubit[], Qubit[]) => Unit is Adj + Ctl |
| 166 | ); |
| 167 | |
| 168 | /// # Summary |
| 169 | /// Performs quantum circuit for one-electron operator (p. 54, eq. 70) |
| 170 | internal function MakeOneElectronOperator( |
| 171 | eigenValues : Double[], |
| 172 | eigenVectors : Double[][], |
| 173 | constants : DoubleFactorizedChemistryConstants |
| 174 | ) : OneElectronOperator { |
| 175 | let data = Mapped(eigenValue -> [eigenValue < 0.0], eigenValues); |
| 176 | let prepare = MakePrepareArbitrarySuperpositionWithData(constants::TargetError, eigenValues, data); |
| 177 | |
| 178 | OneElectronOperator( |
| 179 | prepare::NGarbageQubits, |
| 180 | OneElectronOperatorOperation(eigenVectors, constants, prepare, _, _, _, _, _) |
| 181 | ) |
| 182 | } |
| 183 | |
| 184 | internal operation OneElectronOperatorOperation( |
| 185 | eigenVectors : Double[][], |
| 186 | constants : DoubleFactorizedChemistryConstants, |
| 187 | prepare : PrepareArbitrarySuperposition, |
| 188 | register0 : Qubit[], |
| 189 | register1 : Qubit[], |
| 190 | phaseGradientRegister : Qubit[], |
| 191 | offsets : Qubit[], |
| 192 | helper : Qubit[] |
| 193 | ) : Unit is Adj + Ctl { |
| 194 | if BeginEstimateCaching("OneElectronOperator", IsEmpty(offsets) ? 0 | 1) { |
| 195 | // assertions |
| 196 | Fact(Length(helper) == prepare::NGarbageQubits, "invalid number of helper qubits in OneElectronOperatorOperation"); |
| 197 | |
| 198 | let precision = constants::RotationAngleBitPrecision; |
| 199 | let bitstrings = AllEigenVectorsAsBitString(eigenVectors, precision); |
| 200 | |
| 201 | use prepQubits = Qubit[prepare::NIndexQubits]; |
| 202 | use rotationQubits = Qubit[Length(Head(bitstrings))]; |
| 203 | use sign = Qubit(); |
| 204 | use spin = Qubit(); |
| 205 | |
| 206 | within { |
| 207 | if not IsEmpty(offsets) { |
| 208 | prepare::PrepareWithSelect(SelectWithOffset(_, offsets, _, _), prepQubits, [sign], helper); |
| 209 | } else { |
| 210 | prepare::Prepare(prepQubits, [sign], helper); |
| 211 | } |
| 212 | H(spin); |
| 213 | for i in IndexRange(register0) { |
| 214 | Controlled SWAP([spin], (register0[i], register1[i])); |
| 215 | } |
| 216 | if not IsEmpty(offsets) { |
| 217 | RippleCarryCGIncByLE(offsets, prepQubits); |
| 218 | } |
| 219 | Select(bitstrings, prepQubits, rotationQubits); |
| 220 | } apply { |
| 221 | ApplyGivensRotations(phaseGradientRegister, Chunks(precision, rotationQubits), register0); |
| 222 | Z(sign); |
| 223 | } |
| 224 | |
| 225 | EndEstimateCaching(); |
| 226 | } |
| 227 | } |
| 228 | |
| 229 | internal operation SelectWithOffset(data : Bool[][], offset : Qubit[], address : Qubit[], target : Qubit[]) : Unit is Adj + Ctl { |
| 230 | within { |
| 231 | RippleCarryCGIncByLE(offset, address); |
| 232 | } apply { |
| 233 | Select(data, address, target); |
| 234 | } |
| 235 | } |
| 236 | |
| 237 | internal newtype TwoElectronOperator = ( |
| 238 | NGarbageQubits: Int, |
| 239 | Apply: (Qubit[], Qubit[], Qubit[], Qubit[]) => Unit is Ctl |
| 240 | ); |
| 241 | |
| 242 | /// # Summary |
| 243 | /// Performs quantum circuit for two-electron operator (p. 57, eq. 78) |
| 244 | internal function MakeTwoElectronOperator(problem : DoubleFactorizedChemistryProblem, constants : DoubleFactorizedChemistryConstants) : TwoElectronOperator { |
| 245 | let lambdaPrepare = MakePrepareArbitrarySuperposition( |
| 246 | constants::TargetError, |
| 247 | problem::Lambdas |
| 248 | ); |
| 249 | |
| 250 | let eigenValuesFlattened = Flattened(problem::TwoBodyEigenValues); |
| 251 | let eigenVectorsFlattened = Flattened(problem::TwoBodyEigenVectors); |
| 252 | let oneElectronOperator = MakeOneElectronOperator(eigenValuesFlattened, eigenVectorsFlattened, constants); |
| 253 | |
| 254 | let numGarbageQubits = lambdaPrepare::NGarbageQubits + oneElectronOperator::NGarbageQubits; |
| 255 | |
| 256 | TwoElectronOperator(numGarbageQubits, TwoElectronOperatorOperation(problem, lambdaPrepare, oneElectronOperator, _, _, _, _)) |
| 257 | } |
| 258 | |
| 259 | internal operation TwoElectronOperatorOperation( |
| 260 | problem : DoubleFactorizedChemistryProblem, |
| 261 | lambdaPrepare : PrepareArbitrarySuperposition, |
| 262 | oneElectronOperator : OneElectronOperator, |
| 263 | register0 : Qubit[], |
| 264 | register1 : Qubit[], |
| 265 | phaseGradientRegister : Qubit[], |
| 266 | helper : Qubit[] |
| 267 | ) : Unit is Ctl { |
| 268 | let helperParts = Partitioned([lambdaPrepare::NGarbageQubits], helper); |
| 269 | let dataOffsets = ComputeOffsetDataSet(Mapped(Length, problem::TwoBodyEigenValues)); |
| 270 | |
| 271 | use rankQubits = Qubit[lambdaPrepare::NIndexQubits]; |
| 272 | use offsetQubits = Qubit[Length(dataOffsets[0])]; |
| 273 | |
| 274 | within { |
| 275 | lambdaPrepare::Prepare(rankQubits, [], helperParts[0]); |
| 276 | Select(dataOffsets, rankQubits, offsetQubits); |
| 277 | oneElectronOperator::Apply(register0, register1, phaseGradientRegister, offsetQubits, helperParts[1]); |
| 278 | } apply { |
| 279 | ReflectAboutInteger(0, helperParts[1]); |
| 280 | } |
| 281 | } |
| 282 | |
| 283 | internal function ComputeOffsetDataSet(lengths : Int[]) : Bool[][] { |
| 284 | mutable currentOffset = 0; |
| 285 | mutable offsets = [0, size = Length(lengths)]; |
| 286 | |
| 287 | for i in IndexRange(lengths) { |
| 288 | set offsets w/= i <- currentOffset; |
| 289 | set currentOffset += lengths[i]; |
| 290 | } |
| 291 | |
| 292 | let offsetWidth = Ceiling(Lg(IntAsDouble(currentOffset))); |
| 293 | return Mapped(IntAsBoolArray(_, offsetWidth), offsets); |
| 294 | } |
| 295 | |
| 296 | /// # Summary |
| 297 | /// Performs quantum circuit for phase rotations (p. 54, eq. 73) |
| 298 | internal operation ApplyGivensRotations(pgr : Qubit[], rotationQubits : Qubit[][], target : Qubit[]) : Unit is Adj + Ctl { |
| 299 | within { |
| 300 | let windows = Windows(2, target); |
| 301 | for i in RangeReverse(IndexRange(rotationQubits)) { |
| 302 | let rotationControls = rotationQubits[i]; |
| 303 | let qs = windows[i]; |
| 304 | ApplyWithMajoranaCliffords(1, ApplyRotationUsingRippleCarryAddition(pgr, rotationControls, _), qs); |
| 305 | ApplyWithMajoranaCliffords(0, Adjoint ApplyRotationUsingRippleCarryAddition(pgr, rotationControls, _), qs); |
| 306 | } |
| 307 | } apply { |
| 308 | X(Head(target)); |
| 309 | Y(Head(target)); |
| 310 | } |
| 311 | } |
| 312 | |
| 313 | internal operation ApplyRotationUsingRippleCarryAddition(pgr : Qubit[], rotationQubits : Qubit[], target : Qubit) : Unit is Adj { |
| 314 | Controlled Adjoint RippleCarryCGIncByLE([target], (Reversed(Rest(rotationQubits)), pgr)); |
| 315 | } |
| 316 | |
| 317 | internal operation ApplyWithMajoranaCliffords(x : Int, op : (Qubit => Unit is Adj), qs : Qubit[]) : Unit is Adj { |
| 318 | // x is the Majorana type (which can be either 0 or 1, see p. 52, eq. 59) |
| 319 | within { |
| 320 | Adjoint S(qs[x]); |
| 321 | H(qs[x]); |
| 322 | H(qs[1 - x]); |
| 323 | CNOT(qs[1], qs[0]); |
| 324 | } apply { |
| 325 | op(qs[0]); |
| 326 | } |
| 327 | } |
| 328 | |
| 329 | internal function AllEigenVectorsAsBitString(eigenVectors : Double[][], precision : Int) : Bool[][] { |
| 330 | mutable bitstrings = []; |
| 331 | |
| 332 | let tau = 2.0 * PI(); |
| 333 | let preFactor = 2.0^IntAsDouble(precision); |
| 334 | |
| 335 | for eigenVector in eigenVectors { |
| 336 | // Computes rotation angles for Majorana operator ($\vec u$ in p. 52, eq. 55) |
| 337 | mutable result = []; |
| 338 | mutable sins = 1.0; |
| 339 | |
| 340 | for index in 0..Length(eigenVector)-2 { |
| 341 | // We apply MinD, such that rounding errors do not lead to |
| 342 | // an argument for ArcCos which is larger than 1.0. (p. 52, eq. 56) |
| 343 | let theta = sins == 0.0 ? 0.0 | 0.5 * ArcCos(MinD(eigenVector[index] / sins, 1.0)); |
| 344 | |
| 345 | // all angles as bit string |
| 346 | let factor = theta / tau; |
| 347 | set result += Reversed(IntAsBoolArray(IsNaN(factor) ? 0 | Floor(preFactor * factor), precision)); |
| 348 | |
| 349 | set sins *= Sin(2.0 * theta); |
| 350 | } |
| 351 | |
| 352 | set bitstrings += [result]; |
| 353 | } |
| 354 | |
| 355 | bitstrings |
| 356 | } |
| 357 | |
| 358 | internal function IsNaN(value : Double) : Bool { |
| 359 | value != value |
| 360 | } |
| 361 | } |
| 362 | |