Design of higher order response operators

[stigrj]

From @robertodr in #161 :

Some comments and suggestions, mostly due to the fact that I am not quite that fluent in MRChem. The D1/D2 variants are there for ground state and linear response, right? That would mean that for e.g. quadratic you'd have to add a D3 variant, isn't it? Using variadic templates might help handle the code duplication a bit:

template <typename Head, typename... Tail>
class Collection final {
   public:
     Collection(const Head & h, const Tail&... t) : collection_{h, t...} {}
    void print() const {
        for (const auto & f: collection_) {
            f.print();
        }
    }
    private:
    std::array<Head, sizeof...(Tail)+1> collection_;
};

where the type of Head and Tail should be Orbital * in your case. Based on the sizeof...(Tail)+1 (a compile-time constant) you can dispatch all the needed algorithms correctly with not too much hassle (I think!) This is a slightly more verbose solution upon instantiation: you need to give the types of all arguments up front: auto coll = Collection<A, A, A, A>(a1, a2, a3, a4); There's probably a smarter way to coerce the parameter pack to be of one single type though.