Method getK3Local

dscribe/ext/mbtr.cpp:670–825

}

vector<map<string, vector<double>>> MBTR::getK3Local(const vector<int> &indices, const vector<int> &Z, const vector<vector<double>> &distances, const vector<vector<int>> &neighbours, const string &geomFunc, const string &weightFunc, const map<string, double> &parameters, double min, double max, double sigma, int n)
{
    // Initialize some variables outside the loop
    vector<map<string, vector<double>>> k3Maps(indices.size());
    double dx = (max-min)/(n-1);
    double sigmasqrt2 = sigma*sqrt(2.0);
    double start = min-dx/2;
    int iLoc = 0;

    // We loop over the specified indices
    for (const int &i : indices) {
        map<string, vector<double> > k3Map;

        // For each atom we loop only over the atoms triplets that are
        // within the neighbourhood
        const vector<int> &i_neighbours = neighbours[i];
        for (const int &j : i_neighbours) {
            for (const int &k : i_neighbours) {
                // Calculate angle for all index permutations from choosing
                // three out of nAtoms. The same atom cannot be present twice
                // in the permutation.
                if (j != i && k != j && k != i) {

                    // Calculate geometry value
                    double geom;
                    if (geomFunc == "cosine") {
                        geom = k3GeomCosine(i, j, k, distances);
                    } else if (geomFunc == "angle") {
                        geom = k3GeomAngle(i, j, k, distances);
                    } else {
                        throw invalid_argument("Invalid geometry function.");
                    }

                    // Calculate weight value
                    double weight;
                    if (weightFunc == "exp") {
                        double scale = parameters.at("scale");
                        double threshold = parameters.at("threshold");
                        weight = k3WeightExponential(i, j, k, distances, scale);
                        if (weight < threshold) {
                            continue;
                        }
                    } else if (weightFunc == "unity") {
                        weight = k3WeightUnity(i, j, k, distances);
                    } else {
                        throw invalid_argument("Invalid weighting function.");
                    }

                    // Calculate gaussian
                    vector<double> gauss = gaussian(geom, weight, start, dx, sigmasqrt2, n);

                    // Get the index of the present elements in the final vector
                    int i_elem = 0;
                    int j_elem = Z[j];
                    int k_elem = Z[k];
                    int i_index = this->atomicNumberToIndexMap.at(i_elem);
                    int j_index = this->atomicNumberToIndexMap.at(j_elem);
                    int k_index = this->atomicNumberToIndexMap.at(k_elem);