WeightedGrid.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
 
#include "Mathematical/Grid/CubicGrid.hpp"
#include "Mathematical/Grid/Field.hpp"
#include "Mathematical/Representation/Array.hpp"
#include "Mathematical/Representation/Matrix.hpp"
 
 
namespace GQCP {
 
 
class WeightedGrid {
private:
    ArrayX<double> m_weights;  // a 1-D array containing the weights for each of the grid points
 
    std::vector<Vector<double, 3>> m_points;  // the grid points
 
 
public:
    // CONSTRUCTORS
 
    WeightedGrid(const std::vector<Vector<double, 3>>& points, const ArrayX<double>& weights);
 
 
    // NAMED CONSTRUCTORS
 
    static WeightedGrid FromCubicGrid(const CubicGrid& cubic_grid);
 
    static WeightedGrid ReadIntegrationGridFile(const std::string& filename);
 
 
    // PUBLIC METHODS
 
    template <typename T>
    T integrate(const Field<T>& field) const {
 
        // A KISS-implementation of integration: multiplying every field value by the weight of the associated grid point.
        auto result = field.value(0) * this->weight(0);  // this makes sure that 'result' is already initialized to the correct type, e.g. when a Vector is initialized, it doesn't automatically have the required size
        for (size_t i = 1; i < this->size(); i++) {
            result += field.value(i) * this->weight(i);
        }
 
        return result;
    }
 
 
    size_t numberOfPoints() const { return this->m_points.size(); }
 
    const Vector<double, 3>& point(const size_t index) const { return this->m_points[index]; }
 
    const std::vector<Vector<double, 3>>& points() const { return this->m_points; }
 
    size_t size() const { return this->numberOfPoints(); }
 
    double weight(const size_t index) const { return this->m_weights(index); }
 
    const ArrayX<double>& weights() const { return this->m_weights; };
};
 
 
}  // namespace GQCP