# Custom scalar types and parametric_loperator

The goal of this example is twofold. First, it shows how to use user-defined types as ScalarType’s, i.e. types of coefficient in expressions. Second, it demonstrates functionality of the parametric_loperator class.

As for the custom scalar type, we will define a class implementing the algebra of polynomials

$P_n(x) = C_0 + C_1 x + C_2 x^2 + \ldots + C_n x^n.$

$$x$$ can be understood as a parameter (for instance, time or an external field) a quantum-mechanical operator depends on. One could envision other practically relevant choices of ScalarType such as interpolators based on the Chebyshev polynomial expansion or Padé approximants. The expression constructed in the program is the Hamiltonian of a displaced quantum harmonic oscillator with frequency $$\omega_0 = 2$$ and the displacement depending on an external parameter $$\lambda$$, $$g(\lambda) = 3\lambda$$.

$\hat H = \frac{\omega_0}{2} + \omega_0 \left(a^\dagger - g(\lambda)\right) \left(a - g(\lambda)\right).$

An instance of the parametric_loperator class is constructed out of this expression and applied to the vacuum state for a few values of $$\lambda$$.

$\hat H|0\rangle = \omega_0\left(g^2(\lambda) + \frac{1}{2}\right)|0\rangle - \omega_0 g(\lambda) |1\rangle.$

The crucial point here is that the slow task of constructing parametric_loperator is performed only once. Substitution of $$\lambda$$ values is delayed until the moment the operator acts on the state.

  1#include <algorithm>
2#include <initializer_list>
3#include <iostream>
4#include <vector>
5
6#include <libcommute/libcommute.hpp>
7
8using namespace libcommute;
9
10// Polynomial of one real variable - a (very) basic implementation.
11// For production code, it is strongly advised to use optimized libraries
12// such as 'polynomial' from Boost Math Toolkit.
13class polynomial {
14
15  // Polynomial coefficients {C_0, C_1, C_2, ..., C_n}.
16  // P_n(x) = C_0 + C_1 x + C_2 x^2 + ... + C_n * x^n.
17  std::vector<double> coefficients;
18
19public:
20  // Construct zero polynomial of degree n
21  explicit polynomial(std::size_t n) : coefficients(n + 1, 0) {}
22
23  // Construct from a list of coefficients
24  explicit polynomial(std::initializer_list<double> coeffs)
25    : coefficients(coeffs) {}
26
27  // Degree of this polynomial
28  std::size_t degree() const { return coefficients.size() - 1; }
29
31  std::vector<double> const& coeffs() const { return coefficients; }
32
33  // Evaluate P_n(x)
34  double operator()(double x) const {
35    double sum = 0;
36    double x_n = 1;
37    for(double c : coefficients) {
38      sum += c * x_n;
39      x_n *= x;
40    }
41    return sum;
42  }
43
45  polynomial operator+(polynomial const& p) const {
46    auto max_degree = std::max(degree(), p.degree());
47    polynomial result(max_degree);
48    for(std::size_t i = 0; i <= max_degree; ++i) {
49      result.coefficients[i] = (i <= degree() ? coefficients[i] : 0) +
50                               (i <= p.degree() ? p.coefficients[i] : 0);
51    }
52    return result;
53  }
54
55  // Subtraction of polynomials
56  polynomial operator-(polynomial const& p) const {
57    auto max_degree = std::max(degree(), p.degree());
58    polynomial result(max_degree);
59    for(std::size_t i = 0; i <= max_degree; ++i) {
60      result.coefficients[i] = (i <= degree() ? coefficients[i] : 0) -
61                               (i <= p.degree() ? p.coefficients[i] : 0);
62    }
63    return result;
64  }
65
66  // Multiplication of polynomials
67  polynomial operator*(polynomial const& p) const {
68    polynomial result(degree() + p.degree());
69    for(std::size_t i = 0; i <= degree(); ++i) {
70      for(std::size_t j = 0; j <= p.degree(); ++j) {
71        result.coefficients[i + j] += coefficients[i] * p.coefficients[j];
72      }
73    }
74    return result;
75  }
76
77  // Multiplication by a constant
78  polynomial operator*(double x) const {
79    polynomial res(*this);
80    std::for_each(res.coefficients.begin(),
81                  res.coefficients.end(),
82                  [x](double& c) { c *= x; });
83    return res;
84  }
85
86  // Multiplication by a constant pre-factor
87  friend polynomial operator*(double x, polynomial const& p) { return p * x; }
88};
89
90// Specialize struct scalar_traits to let libcommute know about our polynomial
91// type so that it can be used as a coefficient type of libcommute::expression.
92namespace libcommute {
93
94template <> struct scalar_traits<polynomial> {
95  // Zero value test
96  static bool is_zero(polynomial const& p) {
97    return std::all_of(p.coeffs().begin(), p.coeffs().end(), [](double c) {
98      return c == 0;
99    });
100  }
101  // Make a constant polynomial from a double value
102  static polynomial make_const(double x) { return polynomial{x}; }
103};
104
105} // namespace libcommute
106
107int main() {
108
109  //
110  // Displaced harmonic oscillator
111  //
112
113  // Oscillator frequency \omega_0 = 2.
114  polynomial w0{2};
115  // Displacement as a linear function of an external parameter,
116  // g(\lambda) = 3\lambda.
117  polynomial g{0, 3};
118
119  // Hamiltonian of the oscillator
120  using namespace static_indices; // For a_dag() and a()
121  auto H = 0.5 * w0 + w0 * (a_dag<polynomial>() - g) * (a<polynomial>() - g);
122
123  // H will act in the Hilbert space 'hs'. Since we are working with a boson,
124  // we have to truncate the Hilbert space and allow a finite number of
125  // excitations (here, N = 0, 1, ..., 2^4 - 1).
126  auto hs = make_hilbert_space(H, boson_es_constructor(4));
127
128  // Parametric quantum operator (depends on parameter \lambda)
129  auto Hop = make_param_loperator(H, hs);
130
131  //
132  // Check that
133  // H|0> = \omega_0 (1/2 + g(\lambda)^2) |0> - \omega_0 g(\lambda) |1>
134  // for a few values of \lambda.
135  //
136
137  // Initial state |0>
138  std::vector<double> ket0(hs.dim());
139  ket0[0] = 1.0;
140
141  for(double lambda : {.0, 1.0, 2.0, 3.0}) {
142    // Act with H(\lambda) upon |0>
143    auto ket = Hop(ket0, lambda);
144
145    // Print all non-zero elements of |ket>
146    std::cout << "\\lambda = " << lambda << " :|ket> = ";
147    for(unsigned int i = 0; i < ket.size(); ++i) {
148      if(ket[i] != 0) std::cout << " +(" << ket[i] << ")|" << i << ">";
149    }
150    std::cout << '\n';
151  }
152
153  return 0;
154}