Method getK2Local

dscribe/ext/mbtr.cpp:584–668

}

vector<map<string, vector<double>>> MBTR::getK2Local(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
    int nPos = indices.size();
    vector<map<string, vector<double> > > k2Maps(nPos);
    double dx = (max-min)/(n-1);
    double sigmasqrt2 = sigma*sqrt(2.0);
    double start = min-dx/2;

    // We loop over the specified indices
    for (int i = 0; i < nPos; ++i) {
        int iTrue = indices[i];
        map<string, vector<double> > k2Map;

        // For each atom we loop only over the neighbours
        const vector<int> &i_neighbours = neighbours[i];
        for (const int &j : i_neighbours) {

            // Self-distances are not considered
            if (iTrue == j) {
                continue;
            }

            // Calculate geometry value
            double geom;
            if (geomFunc == "inverse_distance") {
                geom = k2GeomInverseDistance(i, j, distances);
            } else if (geomFunc == "distance") {
                geom = k2GeomDistance(i, j, 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 = k2WeightExponential(i, j, distances, scale);
                if (weight < threshold) {
                    continue;
                }
            } else if (weightFunc == "unity") {
                weight = k2WeightUnity(i, j, 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 i_index = this->atomicNumberToIndexMap.at(i_elem);
            int j_index = this->atomicNumberToIndexMap.at(j_elem);

            // Save information in the part where j_index >= i_index