colvarcomp_gpath.cpp

// -*- c++ -*-

// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/Colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.

#include <numeric>
#include <algorithm>
#include <cmath>
#include <limits>
#include <fstream>

#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvar.h"
#include "colvarcomp.h"



colvar::CartesianBasedPath::CartesianBasedPath()
{
    x.type(colvarvalue::type_scalar);
    // Don't use implicit gradient
    enable(f_cvc_explicit_gradient);
}


int colvar::CartesianBasedPath::init(std::string const &conf)
{
    int error_code = cvc::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    // Parse selected atoms
    atoms = parse_group(conf, "atoms");
    has_user_defined_fitting = false;
    std::string fitting_conf;
    if (key_lookup(conf, "fittingAtoms", &fitting_conf)) {
        has_user_defined_fitting = true;
    }
    // Lookup reference column of PDB
    // Copied from the RMSD class
    std::string reference_column;
    double reference_column_value = 0.0;
    if (get_keyval(conf, "refPositionsCol", reference_column, std::string(""))) {
        bool found = get_keyval(conf, "refPositionsColValue", reference_column_value, reference_column_value);
        if (found && reference_column_value == 0.0) {
            return cvm::error("Error: refPositionsColValue, if provided, must be non-zero.\n",
                              COLVARS_INPUT_ERROR);
        }
    }
    // Lookup all reference frames
    bool has_frames = true;
    total_reference_frames = 0;
    while (has_frames) {
        std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(total_reference_frames + 1);
        if (key_lookup(conf, reference_position_file_lookup.c_str())) {
            std::string reference_position_filename;
            get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
            std::vector<cvm::atom_pos> reference_position(atoms->size());
            cvm::load_coords(reference_position_filename.c_str(), &reference_position, atoms, reference_column, reference_column_value);
            reference_frames.push_back(reference_position);
            ++total_reference_frames;
        } else {
            has_frames = false;
        }
    }
    // Setup alignment to compute RMSD with respect to reference frames
    for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
        cvm::atom_group* tmp_atoms = parse_group(conf, "atoms");
        if (!has_user_defined_fitting) {
            // Swipe from the rmsd class
            tmp_atoms->enable(f_ag_center);
            tmp_atoms->enable(f_ag_rotate);
            tmp_atoms->ref_pos = reference_frames[i_frame];
            tmp_atoms->center_ref_pos();
            tmp_atoms->enable(f_ag_fit_gradients);
        } else {
            // parse a group of atoms for fitting
            std::string fitting_group_name = std::string("fittingAtoms") + cvm::to_str(i_frame);
            cvm::atom_group* tmp_fitting_atoms = new cvm::atom_group(fitting_group_name.c_str());
            tmp_fitting_atoms->parse(fitting_conf);
            tmp_fitting_atoms->disable(f_ag_scalable);
            tmp_fitting_atoms->fit_gradients.assign(tmp_fitting_atoms->size(), cvm::atom_pos(0.0, 0.0, 0.0));
            std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(i_frame + 1);
            std::string reference_position_filename;
            get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
            std::vector<cvm::atom_pos> reference_fitting_position(tmp_fitting_atoms->size());
            cvm::load_coords(reference_position_filename.c_str(), &reference_fitting_position, tmp_fitting_atoms, reference_column, reference_column_value);
            // setup the atom group for calculating
            tmp_atoms->enable(f_ag_center);
            tmp_atoms->enable(f_ag_rotate);
            tmp_atoms->b_user_defined_fit = true;
            tmp_atoms->disable(f_ag_scalable);
            tmp_atoms->ref_pos = reference_fitting_position;
            tmp_atoms->center_ref_pos();
            tmp_atoms->enable(f_ag_fit_gradients);
            tmp_atoms->enable(f_ag_fitting_group);
            tmp_atoms->fitting_group = tmp_fitting_atoms;
            reference_fitting_frames.push_back(reference_fitting_position);
        }
        tmp_atoms->setup_rotation_derivative();
        comp_atoms.push_back(tmp_atoms);
    }

    return error_code;
}

colvar::CartesianBasedPath::~CartesianBasedPath() {
    for (auto it_comp_atoms = comp_atoms.begin(); it_comp_atoms != comp_atoms.end(); ++it_comp_atoms) {
        if (*it_comp_atoms != nullptr) {
            delete (*it_comp_atoms);
            (*it_comp_atoms) = nullptr;
        }
    }
    // Avoid double-freeing due to CVC-in-CVC construct
    atom_groups.clear();
}

void colvar::CartesianBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& result) {
    for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
        cvm::real frame_rmsd = 0.0;
        for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            frame_rmsd += ((*(comp_atoms[i_frame]))[i_atom].pos - reference_frames[i_frame][i_atom]).norm2();
        }
        frame_rmsd /= cvm::real(atoms->size());
        frame_rmsd = cvm::sqrt(frame_rmsd);
        result[i_frame] = frame_rmsd;
    }
}

// mainly used for determining the lambda value for arithmetic path
void colvar::CartesianBasedPath::computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& result) {
    for (size_t i_frame = 0; i_frame < reference_frames.size() - 1; ++i_frame) {
        std::vector<cvm::atom_pos> this_frame_atom_pos(reference_frames[i_frame].size());
        std::vector<cvm::atom_pos> next_frame_atom_pos(reference_frames[i_frame + 1].size());
        cvm::real frame_rmsd = 0.0;
        const size_t this_index = i_frame;
        const size_t next_index = i_frame + 1;
        // compute COM of two successive images, respectively
        cvm::atom_pos reference_cog_this, reference_cog_next;
        for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            reference_cog_this += reference_frames[this_index][i_atom];
            reference_cog_next += reference_frames[next_index][i_atom];
        }
        reference_cog_this /= reference_frames[this_index].size();
        reference_cog_next /= reference_frames[next_index].size();
        // move all atoms to COM
        for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            this_frame_atom_pos[i_atom] = reference_frames[this_index][i_atom] - reference_cog_this;
            next_frame_atom_pos[i_atom] = reference_frames[next_index][i_atom] - reference_cog_next;
        }
        cvm::rotation rot_this_to_next;
        // compute the optimal rotation
        rot_this_to_next.calc_optimal_rotation(this_frame_atom_pos, next_frame_atom_pos);
        // compute rmsd between reference frames
        for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            frame_rmsd += (rot_this_to_next.q.rotate(this_frame_atom_pos[i_atom]) - next_frame_atom_pos[i_atom]).norm2();
        }
        frame_rmsd /= cvm::real(atoms->size());
        frame_rmsd = cvm::sqrt(frame_rmsd);
        result[i_frame] = frame_rmsd;
    }
}


void colvar::CartesianBasedPath::apply_force(colvarvalue const &force)
{
    cvm::error("Error: using apply_force() in a component of type CartesianBasedPath, which is abstract.\n",
               COLVARS_BUG_ERROR);
}



colvar::gspath::gspath()
{
    set_function_type("gspath");
}


int colvar::gspath::init(std::string const &conf)
{
    int error_code = CartesianBasedPath::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
    if (use_second_closest_frame == true) {
        cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
    } else {
        cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
    }
    get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
    if (use_third_closest_frame == true) {
        cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
    } else {
        cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
    }
    if (total_reference_frames < 2) {
        return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) +
                              " reference frames, but gspath requires at least 2 frames.\n",
                          COLVARS_INPUT_ERROR);
    }
    GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::S>::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame);
    cvm::log(std::string("Geometric pathCV(s) is initialized.\n"));
    cvm::log(std::string("Geometric pathCV(s) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
    return error_code;
}

void colvar::gspath::updateDistanceToReferenceFrames() {
    computeDistanceToReferenceFrames(frame_distances);
}

void colvar::gspath::prepareVectors() {
    size_t i_atom;
    for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        // v1 = s_m - z
        v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
        // v2 = z - s_(m-1)
        v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
    }
    if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
        cvm::atom_pos reference_cog_1, reference_cog_2;
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
            reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
        }
        reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
        reference_cog_2 /= cvm::real(reference_frames[min_frame_index_2].size());
        std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
        std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
            tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
        }
        if (has_user_defined_fitting) {
            cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
                reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
                reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
            }
            reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
            reference_fitting_cog_2 /= cvm::real(reference_fitting_frames[min_frame_index_2].size());
            std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
            std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2(reference_fitting_frames[min_frame_index_2].size());
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
                tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
                tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
            }
            rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
        } else {
            rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
        }
        const auto rot_mat_v3 = rot_v3.matrix();
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            v3[i_atom] = rot_mat_v3 * tmp_reference_frame_1[i_atom] - tmp_reference_frame_2[i_atom];
        }
    } else {
        cvm::atom_pos reference_cog_1, reference_cog_3;
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
            reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
        }
        reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
        reference_cog_3 /= cvm::real(reference_frames[min_frame_index_3].size());
        std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
        std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
            tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
        }
        if (has_user_defined_fitting) {
            cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_3;
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
                reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
                reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
            }
            reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
            reference_fitting_cog_3 /= cvm::real(reference_fitting_frames[min_frame_index_3].size());
            std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
            std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
                tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
                tmp_reference_fitting_frame_3[i_atom] = reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
            }
            rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
        } else {
            rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
        }
        const auto rot_mat_v3 = rot_v3.matrix();
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            // v3 = s_(m+1) - s_m
            v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_mat_v3 * tmp_reference_frame_1[i_atom];
        }
    }
}

void colvar::gspath::calc_value() {
    computeValue();
    x = s;
}

void colvar::gspath::calc_gradients() {
    computeDerivatives();
    cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
    // dS(v1, v2(r), v3) / dr = ∂S/∂v1 * dv1/dr + ∂S/∂v2 * dv2/dr
    // dv1/dr = [fitting matrix 1][-1, ..., -1]
    // dv2/dr = [fitting matrix 2][1, ..., 1]
    // ∂S/∂v1 = ± (∂f/∂v1) / (2M)
    // ∂S/∂v2 = ± (∂f/∂v2) / (2M)
    // dS(v1, v2(r), v3) / dr = -1.0 * ± (∂f/∂v1) / (2M) + ± (∂f/∂v2) / (2M)
    for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        tmp_atom_grad_v1[0] = -1.0 * sign * 0.5 * dfdv1[i_atom][0] / M;
        tmp_atom_grad_v1[1] = -1.0 * sign * 0.5 * dfdv1[i_atom][1] / M;
        tmp_atom_grad_v1[2] = -1.0 * sign * 0.5 * dfdv1[i_atom][2] / M;
        tmp_atom_grad_v2[0] = sign * 0.5 * dfdv2[i_atom][0] / M;
        tmp_atom_grad_v2[1] = sign * 0.5 * dfdv2[i_atom][1] / M;
        tmp_atom_grad_v2[2] = sign * 0.5 * dfdv2[i_atom][2] / M;
        (*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
        (*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
    }
}

void colvar::gspath::apply_force(colvarvalue const &force) {
    // The force applied to this CV is scalar type
    cvm::real const &F = force.real_value;
    (*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
    (*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
}


colvar::gzpath::gzpath()
{
    set_function_type("gzpath");
}

int colvar::gzpath::init(std::string const &conf)
{
    int error_code = CartesianBasedPath::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
    if (use_second_closest_frame == true) {
        cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
    } else {
        cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
    }
    get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
    if (use_third_closest_frame == true) {
        cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
    } else {
        cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
    }
    bool b_use_z_square = false;
    get_keyval(conf, "useZsquare", b_use_z_square, false);
    if (b_use_z_square == true) {
        cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
    }
    if (total_reference_frames < 2) {
        return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) +
                               " reference frames, but gzpath requires at least 2 frames.\n",
                          COLVARS_INPUT_ERROR);
    }
    GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::Z>::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square);
    // Logging
    cvm::log(std::string("Geometric pathCV(z) is initialized.\n"));
    cvm::log(std::string("Geometric pathCV(z) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
    return error_code;
}

void colvar::gzpath::updateDistanceToReferenceFrames() {
    computeDistanceToReferenceFrames(frame_distances);
}

void colvar::gzpath::prepareVectors() {
    cvm::atom_pos reference_cog_1, reference_cog_2;
    size_t i_atom;
    for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
        reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
    }
    reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
    reference_cog_2 /= cvm::real(reference_frames[min_frame_index_2].size());
    std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
    std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
    for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
        tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
    }
    std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1;
    std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2;
    if (has_user_defined_fitting) {
        cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
        for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
            reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
            reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
        }
        reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
        reference_fitting_cog_2 /= cvm::real(reference_fitting_frames[min_frame_index_2].size());
        tmp_reference_fitting_frame_1.resize(reference_fitting_frames[min_frame_index_1].size());
        tmp_reference_fitting_frame_2.resize(reference_fitting_frames[min_frame_index_2].size());
        for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
            tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
            tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
        }
        rot_v4.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
    } else {
        rot_v4.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
    }
    const auto rot_mat_v4 = rot_v4.matrix();
    for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
        v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
        // v4 only computes in gzpath
        // v4 = s_m - s_(m-1)
        v4[i_atom] = rot_mat_v4 * tmp_reference_frame_1[i_atom] - tmp_reference_frame_2[i_atom];
    }
    if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
        v3 = v4;
    } else {
        cvm::atom_pos reference_cog_3;
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
        }
        reference_cog_3 /= cvm::real(reference_frames[min_frame_index_3].size());
        std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
        }
        if (has_user_defined_fitting) {
            cvm::atom_pos reference_fitting_cog_3;
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
                reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
            }
            reference_fitting_cog_3 /= cvm::real(reference_fitting_frames[min_frame_index_3].size());
            std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
            for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
                tmp_reference_fitting_frame_3[i_atom] =  reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
            }
            rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
        } else {
            rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
        }
        const auto rot_mat_v3 = rot_v3.matrix();
        for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
            // v3 = s_(m+1) - s_m
            v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_mat_v3 * tmp_reference_frame_1[i_atom];
        }
    }
}

void colvar::gzpath::calc_value() {
    computeValue();
    x = z;
}

void colvar::gzpath::calc_gradients() {
    computeDerivatives();
    cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
    for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
        tmp_atom_grad_v1 = -1.0 * dzdv1[i_atom];
        tmp_atom_grad_v2 = dzdv2[i_atom];
        (*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
        (*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
    }
}

void colvar::gzpath::apply_force(colvarvalue const &force) {
    // The force applied to this CV is scalar type
    cvm::real const &F = force.real_value;
    (*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
    (*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
}


colvar::CVBasedPath::CVBasedPath()
{
    set_function_type("gspathCV");
    x.type(colvarvalue::type_scalar);
}


int colvar::CVBasedPath::init(std::string const &conf)
{
    int error_code = cvc::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    // Lookup all available sub-cvcs
    for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) {
        if (key_lookup(conf, it_cv_map->first.c_str())) {
            std::vector<std::string> sub_cvc_confs;
            get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs);
            for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) {
                cv.push_back((it_cv_map->second)());
                cv.back()->init(*(it_sub_cvc_conf));
            }
        }
    }
    // Sort all sub CVs by their names
    std::sort(cv.begin(), cv.end(), colvar::compare_cvc);
    // Register atom groups and determine the colvar type for reference
    std::vector<colvarvalue> tmp_cv;
    for (auto it_sub_cv = cv.begin(); it_sub_cv != cv.end(); ++it_sub_cv) {
        for (auto it_atom_group = (*it_sub_cv)->atom_groups.begin(); it_atom_group != (*it_sub_cv)->atom_groups.end(); ++it_atom_group) {
            register_atom_group(*it_atom_group);
        }
        colvarvalue tmp_i_cv((*it_sub_cv)->value());
        tmp_i_cv.reset();
        tmp_cv.push_back(tmp_i_cv);
    }
    // Read path file
    // Lookup all reference CV values
    std::string path_filename;
    get_keyval(conf, "pathFile", path_filename);
    cvm::log(std::string("Reading path file: ") + path_filename + std::string("\n"));
    auto &ifs_path = cvm::main()->proxy->input_stream(path_filename);
    if (!ifs_path) {
        return COLVARS_INPUT_ERROR;
    }
    std::string line;
    const std::string token(" ");
    total_reference_frames = 0;
    while (std::getline(ifs_path, line)) {
        std::vector<std::string> fields;
        split_string(line, token, fields);
        size_t num_value_required = 0;
        cvm::log(std::string("Reading reference frame ") + cvm::to_str(total_reference_frames + 1) + std::string("\n"));
        for (size_t i_cv = 0; i_cv < tmp_cv.size(); ++i_cv) {
            const size_t value_size = tmp_cv[i_cv].size();
            num_value_required += value_size;
            cvm::log(std::string("Reading CV ") + cv[i_cv]->name + std::string(" with ") + cvm::to_str(value_size) + std::string(" value(s)\n"));
            if (num_value_required <= fields.size()) {
                size_t start_index = num_value_required - value_size;
                for (size_t i = start_index; i < num_value_required; ++i) {
                    tmp_cv[i_cv][i - start_index] = std::atof(fields[i].c_str());
                    cvm::log(cvm::to_str(tmp_cv[i_cv][i - start_index]));
                }
            } else {
                error_code = cvm::error("Error: incorrect format of path file.\n", COLVARS_INPUT_ERROR);
                return error_code;
            }
        }
        if (!fields.empty()) {
            ref_cv.push_back(tmp_cv);
            ++total_reference_frames;
        }
    }
    cvm::main()->proxy->close_input_stream(path_filename);
    if (total_reference_frames <= 1) {
      error_code = cvm::error(
          "Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n",
          COLVARS_INPUT_ERROR);
      return error_code;
    }
    if (cv.size() == 0) {
      error_code =
          cvm::error("Error: the CV " + name + " expects one or more nesting components.\n",
                     COLVARS_INPUT_ERROR);
      return error_code;
    }

    use_explicit_gradients = true;
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
            use_explicit_gradients = false;
        }
    }
    if (!use_explicit_gradients) {
        disable(f_cvc_explicit_gradient);
    }

    return error_code;
}

void colvar::CVBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& result) {
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        cv[i_cv]->calc_value();
    }
    for (size_t i_frame = 0; i_frame < ref_cv.size(); ++i_frame) {
        cvm::real rmsd_i = 0.0;
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]);
            colvarvalue current_cv_value(cv[i_cv]->value());
            // polynomial combination allowed
            if (current_cv_value.type() == colvarvalue::type_scalar) {
                // wrapping is already in dist2
                rmsd_i += cv[i_cv]->dist2(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value);
            } else {
                rmsd_i += cv[i_cv]->dist2(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
            }
        }
        rmsd_i /= cvm::real(cv.size());
        rmsd_i = cvm::sqrt(rmsd_i);
        result[i_frame] = rmsd_i;
    }
}

void colvar::CVBasedPath::computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& result) const {
    if (ref_cv.size() < 2) return;
    for (size_t i_frame = 1; i_frame < ref_cv.size(); ++i_frame) {
        cvm::real dist_ij = 0.0;
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]);
            colvarvalue prev_ref_cv_value(ref_cv[i_frame-1][i_cv]);
            dist_ij += cv[i_cv]->dist2(ref_cv_value, prev_ref_cv_value);
        }
        dist_ij = cvm::sqrt(dist_ij);
        result[i_frame-1] = dist_ij;
    }
}

cvm::real colvar::CVBasedPath::getPolynomialFactorOfCVGradient(size_t i_cv) const {
    cvm::real factor_polynomial = 1.0;
    if (cv[i_cv]->value().type() == colvarvalue::type_scalar) {
        factor_polynomial = cv[i_cv]->sup_coeff * cv[i_cv]->sup_np * cvm::pow(cv[i_cv]->value().real_value, cv[i_cv]->sup_np - 1);
    } else {
        factor_polynomial = cv[i_cv]->sup_coeff;
    }
    return factor_polynomial;
}

colvar::CVBasedPath::~CVBasedPath() {
    // Recall the steps we initialize the sub-CVCs:
    // 1. Lookup all sub-CVCs and then register the atom groups for sub-CVCs
    //    in their constructors;
    // 2. Iterate over all sub-CVCs, get the pointers of their atom groups
    //    groups, and register again in the parent (current) CVC.
    // That being said, the atom groups become children of the sub-CVCs at
    // first, and then become children of the parent CVC.
    // So, to destruct this class (parent CVC class), we need to remove the
    // dependencies of the atom groups to the parent CVC at first.
    remove_all_children();
    // Then we remove the dependencies of the atom groups to the sub-CVCs
    // in their destructors.
    for (auto it = cv.begin(); it != cv.end(); ++it) {
        delete (*it);
    }
    // The last step is cleaning up the list of atom groups.
    atom_groups.clear();
}


void colvar::CVBasedPath::apply_force(colvarvalue const &force)
{
    cvm::error("Error: using apply_force() in a component of type CVBasedPath, which is abstract.\n",
               COLVARS_BUG_ERROR);
}


cvm::real colvar::CVBasedPath::dist2(colvarvalue const &x1, colvarvalue const &x2) const
{
  return x1.dist2(x2);
}


colvarvalue colvar::CVBasedPath::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const
{
  return x1.dist2_grad(x2);
}


colvarvalue colvar::CVBasedPath::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const
{
  return x2.dist2_grad(x1);
}


void colvar::CVBasedPath::wrap(colvarvalue & /* x_unwrapped */) const {}



colvar::gspathCV::gspathCV()
{
    set_function_type("gspathCV");
    x.type(colvarvalue::type_scalar);
}

int colvar::gspathCV::init(std::string const &conf)
{
    int error_code = CVBasedPath::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
    // Initialize variables for future calculation
    get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
    if (use_second_closest_frame == true) {
        cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
    } else {
        cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
    }
    get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
    if (use_third_closest_frame == true) {
        cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
    } else {
        cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
    }
    if (total_reference_frames < 2) {
        return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspathCV requires at least 2 frames.\n", COLVARS_INPUT_ERROR);
    }
    GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame);
    return error_code;
}

colvar::gspathCV::~gspathCV() {}

void colvar::gspathCV::updateDistanceToReferenceFrames() {
    computeDistanceToReferenceFrames(frame_distances);
}

void colvar::gspathCV::prepareVectors() {
    // Compute v1, v2 and v3
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // values of sub-cvc are computed in update_distances
        // cv[i_cv]->calc_value();
        colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
        colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
        colvarvalue current_cv_value(cv[i_cv]->value());
        // polynomial combination allowed
        if (current_cv_value.type() == colvarvalue::type_scalar) {
            v1[i_cv] = f1_ref_cv_i_value.real_value - cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
            v2[i_cv] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)) - f2_ref_cv_i_value.real_value;
        } else {
            v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
            v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
        }
        cv[i_cv]->wrap(v1[i_cv]);
        cv[i_cv]->wrap(v2[i_cv]);
    }
    if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
            cv[i_cv]->wrap(v3[i_cv]);
        }
    } else {
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
            cv[i_cv]->wrap(v3[i_cv]);
        }
    }
}

void colvar::gspathCV::calc_value() {
    computeValue();
    x = s;
}

void colvar::gspathCV::calc_gradients() {
    computeDerivatives();
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // No matter whether the i-th cv uses implicit gradient, compute it first.
        cv[i_cv]->calc_gradients();
        // If the gradient is not implicit, then add the gradients to its atom groups
        if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
            // Temporary variables storing gradients
            colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
            colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
            // Compute factors for polynomial combinations
            cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
            // Loop over all elements of the corresponding colvar value
            for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
                // ds/dz, z = vector of CVs
                tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
                tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
                // Apply the gradients to the atom groups in i-th cv
                // Loop over all atom groups
                for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
                    // Loop over all atoms in the k-th atom group
                    for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
                        // Chain rule
                        (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * ((*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v1[j_elem] + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v2[j_elem]);
                    }
                }
            }
        }
    }
}

void colvar::gspathCV::apply_force(colvarvalue const &force) {
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // If this CV us explicit gradients, then atomic gradients is already calculated
        // We can apply the force to atom groups directly
        if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
            for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
                (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
            }
        } else {
            // Temporary variables storing gradients
            colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
            colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
            // Compute factors for polynomial combinations
            cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
            for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
                // ds/dz, z = vector of CVs
                tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
                tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
            }
            colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
            cv[i_cv]->apply_force(cv_force);
        }
    }
}

colvar::gzpathCV::gzpathCV()
{
    set_function_type("gzpathCV");
}

int colvar::gzpathCV::init(std::string const &conf) {

    int error_code = CVBasedPath::init(conf);
    if (error_code != COLVARS_OK) return error_code;

    cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
    // Initialize variables for future calculation
    M = cvm::real(total_reference_frames - 1);
    m = 1.0;
    get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
    if (use_second_closest_frame == true) {
        cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
    } else {
        cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
    }
    get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
    if (use_third_closest_frame == true) {
        cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
    } else {
        cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
    }
    bool b_use_z_square = false;
    get_keyval(conf, "useZsquare", b_use_z_square, false);
    if (b_use_z_square == true) {
        cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
    }
    if (total_reference_frames < 2) {
        return cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) +
                                     " reference frames, but gzpathCV requires at least 2 frames.\n",
                                 COLVARS_INPUT_ERROR);
    }
    GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square);

    return error_code;
}

colvar::gzpathCV::~gzpathCV() {
}

void colvar::gzpathCV::updateDistanceToReferenceFrames() {
    computeDistanceToReferenceFrames(frame_distances);
}

void colvar::gzpathCV::prepareVectors() {
    // Compute v1, v2 and v3
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // values of sub-cvc are computed in update_distances
        // cv[i_cv]->calc_value();
        colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
        colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
        colvarvalue current_cv_value(cv[i_cv]->value());
        // polynomial combination allowed
        if (current_cv_value.type() == colvarvalue::type_scalar) {
            v1[i_cv] = f1_ref_cv_i_value.real_value - cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np));
            v2[i_cv] = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)) - f2_ref_cv_i_value.real_value;
        } else {
            v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
            v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
        }
        v4[i_cv] = f1_ref_cv_i_value - f2_ref_cv_i_value;
        cv[i_cv]->wrap(v1[i_cv]);
        cv[i_cv]->wrap(v2[i_cv]);
        cv[i_cv]->wrap(v4[i_cv]);
    }
    if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
            cv[i_cv]->wrap(v3[i_cv]);
        }
    } else {
        for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
            v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
            cv[i_cv]->wrap(v3[i_cv]);
        }
    }
}

void colvar::gzpathCV::calc_value() {
    computeValue();
    x = z;
}

void colvar::gzpathCV::calc_gradients() {
    computeDerivatives();
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // No matter whether the i-th cv uses implicit gradient, compute it first.
        cv[i_cv]->calc_gradients();
        // If the gradient is not implicit, then add the gradients to its atom groups
        if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
            // Temporary variables storing gradients
            colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
            colvarvalue tmp_cv_grad_v2 =  1.0 * dzdv2[i_cv];
            // Compute factors for polynomial combinations
            cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
            for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
                // Apply the gradients to the atom groups in i-th cv
                // Loop over all atom groups
                for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
                    // Loop over all atoms in the k-th atom group
                    for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
                        // Chain rule
                        (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = factor_polynomial * ((*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v1[j_elem] + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad * tmp_cv_grad_v2[j_elem]);
                    }
                }
            }
        }
    }
}

void colvar::gzpathCV::apply_force(colvarvalue const &force) {
    for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
        // If this CV us explicit gradients, then atomic gradients is already calculated
        // We can apply the force to atom groups directly
        if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
            for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
                (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
            }
        } else {
            colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
            colvarvalue tmp_cv_grad_v2 =  1.0 * dzdv2[i_cv];
            // Temporary variables storing gradients
            // Compute factors for polynomial combinations
            cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
            colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
            cv[i_cv]->apply_force(cv_force);
        }
    }
}