Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvarcomp_gpath.C

Go to the documentation of this file.
00001 #if (__cplusplus >= 201103L)
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <numeric>
00011 #include <algorithm>
00012 #include <cmath>
00013 #include <cstdlib>
00014 #include <limits>
00015 #include <fstream>
00016 
00017 #include "colvarmodule.h"
00018 #include "colvarvalue.h"
00019 #include "colvarparse.h"
00020 #include "colvar.h"
00021 #include "colvarcomp.h"
00022 
00023 colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(conf), atoms(nullptr), reference_frames(0) {
00024     // Parse selected atoms
00025     atoms = parse_group(conf, "atoms");
00026     has_user_defined_fitting = false;
00027     std::string fitting_conf;
00028     if (key_lookup(conf, "fittingAtoms", &fitting_conf)) {
00029         has_user_defined_fitting = true;
00030     }
00031     // Lookup reference column of PDB
00032     // Copied from the RMSD class
00033     std::string reference_column;
00034     double reference_column_value;
00035     if (get_keyval(conf, "refPositionsCol", reference_column, std::string(""))) {
00036         bool found = get_keyval(conf, "refPositionsColValue", reference_column_value, 0.0);
00037         if (found && reference_column_value == 0.0) {
00038           cvm::error("Error: refPositionsColValue, "
00039                      "if provided, must be non-zero.\n");
00040           return;
00041         }
00042     }
00043     // Lookup all reference frames
00044     bool has_frames = true;
00045     total_reference_frames = 0;
00046     while (has_frames) {
00047         std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(total_reference_frames + 1);
00048         if (key_lookup(conf, reference_position_file_lookup.c_str())) {
00049             std::string reference_position_filename;
00050             get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
00051             std::vector<cvm::atom_pos> reference_position(atoms->size());
00052             cvm::load_coords(reference_position_filename.c_str(), &reference_position, atoms, reference_column, reference_column_value);
00053             reference_frames.push_back(reference_position);
00054             ++total_reference_frames;
00055         } else {
00056             has_frames = false;
00057         }
00058     }
00059     // Setup alignment to compute RMSD with respect to reference frames
00060     for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
00061         cvm::atom_group* tmp_atoms = parse_group(conf, "atoms");
00062         if (!has_user_defined_fitting) {
00063             // Swipe from the rmsd class
00064             tmp_atoms->enable(f_ag_center);
00065             tmp_atoms->enable(f_ag_rotate);
00066             tmp_atoms->ref_pos = reference_frames[i_frame];
00067             tmp_atoms->center_ref_pos();
00068             tmp_atoms->enable(f_ag_fit_gradients);
00069             tmp_atoms->rot.request_group1_gradients(tmp_atoms->size());
00070             tmp_atoms->rot.request_group2_gradients(tmp_atoms->size());
00071             comp_atoms.push_back(tmp_atoms);
00072         } else {
00073             // parse a group of atoms for fitting
00074             std::string fitting_group_name = std::string("fittingAtoms") + cvm::to_str(i_frame);
00075             cvm::atom_group* tmp_fitting_atoms = new cvm::atom_group(fitting_group_name.c_str());
00076             tmp_fitting_atoms->parse(fitting_conf);
00077             tmp_fitting_atoms->disable(f_ag_scalable);
00078             tmp_fitting_atoms->fit_gradients.assign(tmp_fitting_atoms->size(), cvm::atom_pos(0.0, 0.0, 0.0));
00079             std::string reference_position_file_lookup = "refPositionsFile" + cvm::to_str(i_frame + 1);
00080             std::string reference_position_filename;
00081             get_keyval(conf, reference_position_file_lookup.c_str(), reference_position_filename, std::string(""));
00082             std::vector<cvm::atom_pos> reference_fitting_position(tmp_fitting_atoms->size());
00083             cvm::load_coords(reference_position_filename.c_str(), &reference_fitting_position, tmp_fitting_atoms, reference_column, reference_column_value);
00084             // setup the atom group for calculating
00085             tmp_atoms->enable(f_ag_center);
00086             tmp_atoms->enable(f_ag_rotate);
00087             tmp_atoms->b_user_defined_fit = true;
00088             tmp_atoms->disable(f_ag_scalable);
00089             tmp_atoms->ref_pos = reference_fitting_position;
00090             tmp_atoms->center_ref_pos();
00091             tmp_atoms->enable(f_ag_fit_gradients);
00092             tmp_atoms->enable(f_ag_fitting_group);
00093             tmp_atoms->fitting_group = tmp_fitting_atoms;
00094             tmp_atoms->rot.request_group1_gradients(tmp_fitting_atoms->size());
00095             tmp_atoms->rot.request_group2_gradients(tmp_fitting_atoms->size());
00096             reference_fitting_frames.push_back(reference_fitting_position);
00097             comp_atoms.push_back(tmp_atoms);
00098         }
00099     }
00100     x.type(colvarvalue::type_scalar);
00101     // Don't use implicit gradient
00102     enable(f_cvc_explicit_gradient);
00103 }
00104 
00105 colvar::CartesianBasedPath::~CartesianBasedPath() {
00106     for (auto it_comp_atoms = comp_atoms.begin(); it_comp_atoms != comp_atoms.end(); ++it_comp_atoms) {
00107         if (*it_comp_atoms != nullptr) {
00108             delete (*it_comp_atoms);
00109             (*it_comp_atoms) = nullptr;
00110         }
00111     }
00112     // Avoid double-freeing due to CVC-in-CVC construct
00113     atom_groups.clear();
00114 }
00115 
00116 void colvar::CartesianBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& result) {
00117     for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) {
00118         cvm::real frame_rmsd = 0.0;
00119         for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00120             frame_rmsd += ((*(comp_atoms[i_frame]))[i_atom].pos - reference_frames[i_frame][i_atom]).norm2();
00121         }
00122         frame_rmsd /= cvm::real(atoms->size());
00123         frame_rmsd = cvm::sqrt(frame_rmsd);
00124         result[i_frame] = frame_rmsd;
00125     }
00126 }
00127 
00128 colvar::gspath::gspath(std::string const &conf): CartesianBasedPath(conf) {
00129     set_function_type("gspath");
00130     get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
00131     if (use_second_closest_frame == true) {
00132         cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
00133     } else {
00134         cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
00135     }
00136     get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
00137     if (use_third_closest_frame == true) {
00138         cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
00139     } else {
00140         cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
00141     }
00142     if (total_reference_frames < 2) {
00143         cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspath requires at least 2 frames.\n");
00144         return;
00145     }
00146     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);
00147     cvm::log(std::string("Geometric pathCV(s) is initialized.\n"));
00148     cvm::log(std::string("Geometric pathCV(s) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
00149 }
00150 
00151 void colvar::gspath::updateDistanceToReferenceFrames() {
00152     computeDistanceToReferenceFrames(frame_distances);
00153 }
00154 
00155 void colvar::gspath::prepareVectors() {
00156     size_t i_atom;
00157     for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00158         // v1 = s_m - z
00159         v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
00160         // v2 = z - s_(m-1)
00161         v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
00162     }
00163     if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
00164         cvm::atom_pos reference_cog_1, reference_cog_2;
00165         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00166             reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
00167             reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
00168         }
00169         reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
00170         reference_cog_2 /= cvm::real(reference_frames[min_frame_index_2].size());
00171         std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
00172         std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
00173         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00174             tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
00175             tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
00176         }
00177         if (has_user_defined_fitting) {
00178             cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
00179             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00180                 reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
00181                 reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
00182             }
00183             reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
00184             reference_fitting_cog_2 /= cvm::real(reference_fitting_frames[min_frame_index_2].size());
00185             std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
00186             std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2(reference_fitting_frames[min_frame_index_2].size());
00187             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00188                 tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
00189                 tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
00190             }
00191             rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
00192         } else {
00193             rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
00194         }
00195         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00196             v3[i_atom] = rot_v3.q.rotate(tmp_reference_frame_1[i_atom]) - tmp_reference_frame_2[i_atom];
00197         }
00198     } else {
00199         cvm::atom_pos reference_cog_1, reference_cog_3;
00200         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00201             reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
00202             reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
00203         }
00204         reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
00205         reference_cog_3 /= cvm::real(reference_frames[min_frame_index_3].size());
00206         std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
00207         std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
00208         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00209             tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
00210             tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
00211         }
00212         if (has_user_defined_fitting) {
00213             cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_3;
00214             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00215                 reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
00216                 reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
00217             }
00218             reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
00219             reference_fitting_cog_3 /= cvm::real(reference_fitting_frames[min_frame_index_3].size());
00220             std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1(reference_fitting_frames[min_frame_index_1].size());
00221             std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
00222             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00223                 tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
00224                 tmp_reference_fitting_frame_3[i_atom] = reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
00225             }
00226             rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
00227         } else {
00228             rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
00229         }
00230         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00231             // v3 = s_(m+1) - s_m
00232             v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_v3.q.rotate(tmp_reference_frame_1[i_atom]);
00233         }
00234     }
00235 }
00236 
00237 void colvar::gspath::calc_value() {
00238     computeValue();
00239     x = s;
00240 }
00241 
00242 void colvar::gspath::calc_gradients() {
00243     computeDerivatives();
00244     cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
00245     // dS(v1, v2(r), v3) / dr = ∂S/∂v1 * dv1/dr + ∂S/∂v2 * dv2/dr
00246     // dv1/dr = [fitting matrix 1][-1, ..., -1]
00247     // dv2/dr = [fitting matrix 2][1, ..., 1]
00248     // ∂S/∂v1 = ± (∂f/∂v1) / (2M)
00249     // ∂S/∂v2 = ± (∂f/∂v2) / (2M)
00250     // dS(v1, v2(r), v3) / dr = -1.0 * ± (∂f/∂v1) / (2M) + ± (∂f/∂v2) / (2M)
00251     for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00252         tmp_atom_grad_v1[0] = -1.0 * sign * 0.5 * dfdv1[i_atom][0] / M;
00253         tmp_atom_grad_v1[1] = -1.0 * sign * 0.5 * dfdv1[i_atom][1] / M;
00254         tmp_atom_grad_v1[2] = -1.0 * sign * 0.5 * dfdv1[i_atom][2] / M;
00255         tmp_atom_grad_v2[0] = sign * 0.5 * dfdv2[i_atom][0] / M;
00256         tmp_atom_grad_v2[1] = sign * 0.5 * dfdv2[i_atom][1] / M;
00257         tmp_atom_grad_v2[2] = sign * 0.5 * dfdv2[i_atom][2] / M;
00258         (*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
00259         (*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
00260     }
00261 }
00262 
00263 void colvar::gspath::apply_force(colvarvalue const &force) {
00264     // The force applied to this CV is scalar type
00265     cvm::real const &F = force.real_value;
00266     (*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
00267     (*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
00268 }
00269 
00270 colvar::gzpath::gzpath(std::string const &conf): CartesianBasedPath(conf) {
00271     set_function_type("gzpath");
00272     get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
00273     if (use_second_closest_frame == true) {
00274         cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
00275     } else {
00276         cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
00277     }
00278     get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
00279     if (use_third_closest_frame == true) {
00280         cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
00281     } else {
00282         cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
00283     }
00284     bool b_use_z_square = false;
00285     get_keyval(conf, "useZsquare", b_use_z_square, false);
00286     if (b_use_z_square == true) {
00287         cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
00288     }
00289     if (total_reference_frames < 2) {
00290         cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpath requires at least 2 frames.\n");
00291         return;
00292     }
00293     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);
00294     // Logging
00295     cvm::log(std::string("Geometric pathCV(z) is initialized.\n"));
00296     cvm::log(std::string("Geometric pathCV(z) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n"));
00297 }
00298 
00299 void colvar::gzpath::updateDistanceToReferenceFrames() {
00300     computeDistanceToReferenceFrames(frame_distances);
00301 }
00302 
00303 void colvar::gzpath::prepareVectors() {
00304     cvm::atom_pos reference_cog_1, reference_cog_2;
00305     size_t i_atom;
00306     for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00307         reference_cog_1 += reference_frames[min_frame_index_1][i_atom];
00308         reference_cog_2 += reference_frames[min_frame_index_2][i_atom];
00309     }
00310     reference_cog_1 /= cvm::real(reference_frames[min_frame_index_1].size());
00311     reference_cog_2 /= cvm::real(reference_frames[min_frame_index_2].size());
00312     std::vector<cvm::atom_pos> tmp_reference_frame_1(reference_frames[min_frame_index_1].size());
00313     std::vector<cvm::atom_pos> tmp_reference_frame_2(reference_frames[min_frame_index_2].size());
00314     for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00315         tmp_reference_frame_1[i_atom] = reference_frames[min_frame_index_1][i_atom] - reference_cog_1;
00316         tmp_reference_frame_2[i_atom] = reference_frames[min_frame_index_2][i_atom] - reference_cog_2;
00317     }
00318     std::vector<cvm::atom_pos> tmp_reference_fitting_frame_1;
00319     std::vector<cvm::atom_pos> tmp_reference_fitting_frame_2;
00320     if (has_user_defined_fitting) {
00321         cvm::atom_pos reference_fitting_cog_1, reference_fitting_cog_2;
00322         for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00323             reference_fitting_cog_1 += reference_fitting_frames[min_frame_index_1][i_atom];
00324             reference_fitting_cog_2 += reference_fitting_frames[min_frame_index_2][i_atom];
00325         }
00326         reference_fitting_cog_1 /= cvm::real(reference_fitting_frames[min_frame_index_1].size());
00327         reference_fitting_cog_2 /= cvm::real(reference_fitting_frames[min_frame_index_2].size());
00328         tmp_reference_fitting_frame_1.resize(reference_fitting_frames[min_frame_index_1].size());
00329         tmp_reference_fitting_frame_2.resize(reference_fitting_frames[min_frame_index_2].size());
00330         for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_1].size(); ++i_atom) {
00331             tmp_reference_fitting_frame_1[i_atom] = reference_fitting_frames[min_frame_index_1][i_atom] - reference_fitting_cog_1;
00332             tmp_reference_fitting_frame_2[i_atom] = reference_fitting_frames[min_frame_index_2][i_atom] - reference_fitting_cog_2;
00333         }
00334         rot_v4.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_2);
00335     } else {
00336         rot_v4.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_2);
00337     }
00338     for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00339         v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
00340         v2[i_atom] = (*(comp_atoms[min_frame_index_2]))[i_atom].pos - reference_frames[min_frame_index_2][i_atom];
00341         // v4 only computes in gzpath
00342         // v4 = s_m - s_(m-1)
00343         v4[i_atom] = rot_v4.q.rotate(tmp_reference_frame_1[i_atom]) - tmp_reference_frame_2[i_atom];
00344     }
00345     if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
00346         v3 = v4;
00347     } else {
00348         cvm::atom_pos reference_cog_3;
00349         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00350             reference_cog_3 += reference_frames[min_frame_index_3][i_atom];
00351         }
00352         reference_cog_3 /= cvm::real(reference_frames[min_frame_index_3].size());
00353         std::vector<cvm::atom_pos> tmp_reference_frame_3(reference_frames[min_frame_index_3].size());
00354         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00355             tmp_reference_frame_3[i_atom] = reference_frames[min_frame_index_3][i_atom] - reference_cog_3;
00356         }
00357         if (has_user_defined_fitting) {
00358             cvm::atom_pos reference_fitting_cog_3;
00359             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
00360                 reference_fitting_cog_3 += reference_fitting_frames[min_frame_index_3][i_atom];
00361             }
00362             reference_fitting_cog_3 /= cvm::real(reference_fitting_frames[min_frame_index_3].size());
00363             std::vector<cvm::atom_pos> tmp_reference_fitting_frame_3(reference_fitting_frames[min_frame_index_3].size());
00364             for (i_atom = 0; i_atom < reference_fitting_frames[min_frame_index_3].size(); ++i_atom) {
00365                 tmp_reference_fitting_frame_3[i_atom] =  reference_fitting_frames[min_frame_index_3][i_atom] - reference_fitting_cog_3;
00366             }
00367             rot_v3.calc_optimal_rotation(tmp_reference_fitting_frame_1, tmp_reference_fitting_frame_3);
00368         } else {
00369             rot_v3.calc_optimal_rotation(tmp_reference_frame_1, tmp_reference_frame_3);
00370         }
00371         for (i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00372             // v3 = s_(m+1) - s_m
00373             v3[i_atom] = tmp_reference_frame_3[i_atom] - rot_v3.q.rotate(tmp_reference_frame_1[i_atom]);
00374         }
00375     }
00376 }
00377 
00378 void colvar::gzpath::calc_value() {
00379     computeValue();
00380     x = z;
00381 }
00382 
00383 void colvar::gzpath::calc_gradients() {
00384     computeDerivatives();
00385     cvm::rvector tmp_atom_grad_v1, tmp_atom_grad_v2;
00386     for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) {
00387         tmp_atom_grad_v1 = -1.0 * dzdv1[i_atom];
00388         tmp_atom_grad_v2 = dzdv2[i_atom];
00389         (*(comp_atoms[min_frame_index_1]))[i_atom].grad += tmp_atom_grad_v1;
00390         (*(comp_atoms[min_frame_index_2]))[i_atom].grad += tmp_atom_grad_v2;
00391     }
00392 }
00393 
00394 void colvar::gzpath::apply_force(colvarvalue const &force) {
00395     // The force applied to this CV is scalar type
00396     cvm::real const &F = force.real_value;
00397     (*(comp_atoms[min_frame_index_1])).apply_colvar_force(F);
00398     (*(comp_atoms[min_frame_index_2])).apply_colvar_force(F);
00399 }
00400 
00401 
00402 colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) {
00403     // Lookup all available sub-cvcs
00404     for (auto it_cv_map = colvar::get_global_cvc_map().begin(); it_cv_map != colvar::get_global_cvc_map().end(); ++it_cv_map) {
00405         if (key_lookup(conf, it_cv_map->first.c_str())) {
00406             std::vector<std::string> sub_cvc_confs;
00407             get_key_string_multi_value(conf, it_cv_map->first.c_str(), sub_cvc_confs);
00408             for (auto it_sub_cvc_conf = sub_cvc_confs.begin(); it_sub_cvc_conf != sub_cvc_confs.end(); ++it_sub_cvc_conf) {
00409                 cv.push_back((it_cv_map->second)(*(it_sub_cvc_conf)));
00410             }
00411         }
00412     }
00413     // Sort all sub CVs by their names
00414     std::sort(cv.begin(), cv.end(), colvar::compare_cvc);
00415     // Register atom groups and determine the colvar type for reference
00416     std::vector<colvarvalue> tmp_cv;
00417     for (auto it_sub_cv = cv.begin(); it_sub_cv != cv.end(); ++it_sub_cv) {
00418         for (auto it_atom_group = (*it_sub_cv)->atom_groups.begin(); it_atom_group != (*it_sub_cv)->atom_groups.end(); ++it_atom_group) {
00419             register_atom_group(*it_atom_group);
00420         }
00421         colvarvalue tmp_i_cv((*it_sub_cv)->value());
00422         tmp_i_cv.reset();
00423         tmp_cv.push_back(tmp_i_cv);
00424     }
00425     // Read path file
00426     // Lookup all reference CV values
00427     std::string path_filename;
00428     get_keyval(conf, "pathFile", path_filename);
00429     cvm::log(std::string("Reading path file: ") + path_filename + std::string("\n"));
00430     auto &ifs_path = cvm::main()->proxy->input_stream(path_filename);
00431     if (!ifs_path) {
00432         return;
00433     }
00434     std::string line;
00435     const std::string token(" ");
00436     total_reference_frames = 0;
00437     while (std::getline(ifs_path, line)) {
00438         std::vector<std::string> fields;
00439         split_string(line, token, fields);
00440         size_t num_value_required = 0;
00441         cvm::log(std::string("Reading reference frame ") + cvm::to_str(total_reference_frames + 1) + std::string("\n"));
00442         for (size_t i_cv = 0; i_cv < tmp_cv.size(); ++i_cv) {
00443             const size_t value_size = tmp_cv[i_cv].size();
00444             num_value_required += value_size;
00445             cvm::log(std::string("Reading CV ") + cv[i_cv]->name + std::string(" with ") + cvm::to_str(value_size) + std::string(" value(s)\n"));
00446             if (num_value_required <= fields.size()) {
00447                 size_t start_index = num_value_required - value_size;
00448                 for (size_t i = start_index; i < num_value_required; ++i) {
00449                     tmp_cv[i_cv][i - start_index] = std::atof(fields[i].c_str());
00450                     cvm::log(cvm::to_str(tmp_cv[i_cv][i - start_index]));
00451                 }
00452             } else {
00453                 cvm::error("Error: incorrect format of path file.\n");
00454                 return;
00455             }
00456         }
00457         if (!fields.empty()) {
00458             ref_cv.push_back(tmp_cv);
00459             ++total_reference_frames;
00460         }
00461     }
00462     cvm::main()->proxy->close_input_stream(path_filename);
00463     if (total_reference_frames <= 1) {
00464         cvm::error("Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n");
00465         return;
00466     }
00467     if (cv.size() == 0) {
00468         cvm::error("Error: the CV " + name +
00469                    " expects one or more nesting components.\n");
00470         return;
00471     }
00472     x.type(colvarvalue::type_scalar);
00473     use_explicit_gradients = true;
00474     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00475         if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00476             use_explicit_gradients = false;
00477         }
00478     }
00479     if (!use_explicit_gradients) {
00480         disable(f_cvc_explicit_gradient);
00481     }
00482 }
00483 
00484 void colvar::CVBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& result) {
00485     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00486         cv[i_cv]->calc_value();
00487     }
00488     for (size_t i_frame = 0; i_frame < ref_cv.size(); ++i_frame) {
00489         cvm::real rmsd_i = 0.0;
00490         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00491             colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]);
00492             colvarvalue current_cv_value(cv[i_cv]->value());
00493             // polynomial combination allowed
00494             if (current_cv_value.type() == colvarvalue::type_scalar) {
00495                 // wrapping is already in dist2
00496                 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);
00497             } else {
00498                 rmsd_i += cv[i_cv]->dist2(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value);
00499             }
00500         }
00501         rmsd_i /= cvm::real(cv.size());
00502         rmsd_i = cvm::sqrt(rmsd_i);
00503         result[i_frame] = rmsd_i;
00504     }
00505 }
00506 
00507 void colvar::CVBasedPath::computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& result) const {
00508     if (ref_cv.size() < 2) return;
00509     for (size_t i_frame = 1; i_frame < ref_cv.size(); ++i_frame) {
00510         cvm::real dist_ij = 0.0;
00511         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00512             colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]);
00513             colvarvalue prev_ref_cv_value(ref_cv[i_frame-1][i_cv]);
00514             dist_ij += cv[i_cv]->dist2(ref_cv_value, prev_ref_cv_value);
00515         }
00516         dist_ij = cvm::sqrt(dist_ij);
00517         result[i_frame-1] = dist_ij;
00518     }
00519 }
00520 
00521 cvm::real colvar::CVBasedPath::getPolynomialFactorOfCVGradient(size_t i_cv) const {
00522     cvm::real factor_polynomial = 1.0;
00523     if (cv[i_cv]->value().type() == colvarvalue::type_scalar) {
00524         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);
00525     } else {
00526         factor_polynomial = cv[i_cv]->sup_coeff;
00527     }
00528     return factor_polynomial;
00529 }
00530 
00531 colvar::CVBasedPath::~CVBasedPath() {
00532     // Recall the steps we initialize the sub-CVCs:
00533     // 1. Lookup all sub-CVCs and then register the atom groups for sub-CVCs
00534     //    in their constructors;
00535     // 2. Iterate over all sub-CVCs, get the pointers of their atom groups
00536     //    groups, and register again in the parent (current) CVC.
00537     // That being said, the atom groups become children of the sub-CVCs at
00538     // first, and then become children of the parent CVC.
00539     // So, to destruct this class (parent CVC class), we need to remove the
00540     // dependencies of the atom groups to the parent CVC at first.
00541     remove_all_children();
00542     // Then we remove the dependencies of the atom groups to the sub-CVCs
00543     // in their destructors.
00544     for (auto it = cv.begin(); it != cv.end(); ++it) {
00545         delete (*it);
00546     }
00547     // The last step is cleaning up the list of atom groups.
00548     atom_groups.clear();
00549 }
00550 
00551 colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) {
00552     set_function_type("gspathCV");
00553     cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
00554     // Initialize variables for future calculation
00555     get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
00556     if (use_second_closest_frame == true) {
00557         cvm::log(std::string("Geometric path s(σ) will use the second closest frame to compute s_(m-1)\n"));
00558     } else {
00559         cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m-1)\n"));
00560     }
00561     get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
00562     if (use_third_closest_frame == true) {
00563         cvm::log(std::string("Geometric path s(σ) will use the third closest frame to compute s_(m+1)\n"));
00564     } else {
00565         cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n"));
00566     }
00567     if (total_reference_frames < 2) {
00568         cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspathCV requires at least 2 frames.\n");
00569         return;
00570     }
00571     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);
00572     x.type(colvarvalue::type_scalar);
00573 }
00574 
00575 colvar::gspathCV::~gspathCV() {}
00576 
00577 void colvar::gspathCV::updateDistanceToReferenceFrames() {
00578     computeDistanceToReferenceFrames(frame_distances);
00579 }
00580 
00581 void colvar::gspathCV::prepareVectors() {
00582     // Compute v1, v2 and v3
00583     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00584         // values of sub-cvc are computed in update_distances
00585         // cv[i_cv]->calc_value();
00586         colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
00587         colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
00588         colvarvalue current_cv_value(cv[i_cv]->value());
00589         // polynomial combination allowed
00590         if (current_cv_value.type() == colvarvalue::type_scalar) {
00591             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));
00592             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;
00593         } else {
00594             v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
00595             v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
00596         }
00597         cv[i_cv]->wrap(v1[i_cv]);
00598         cv[i_cv]->wrap(v2[i_cv]);
00599     }
00600     if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
00601         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00602             v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
00603             cv[i_cv]->wrap(v3[i_cv]);
00604         }
00605     } else {
00606         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00607             v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
00608             cv[i_cv]->wrap(v3[i_cv]);
00609         }
00610     }
00611 }
00612 
00613 void colvar::gspathCV::calc_value() {
00614     computeValue();
00615     x = s;
00616 }
00617 
00618 void colvar::gspathCV::calc_gradients() {
00619     computeDerivatives();
00620     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00621         // No matter whether the i-th cv uses implicit gradient, compute it first.
00622         cv[i_cv]->calc_gradients();
00623         // If the gradient is not implicit, then add the gradients to its atom groups
00624         if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00625             // Temporary variables storing gradients
00626             colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
00627             colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
00628             // Compute factors for polynomial combinations
00629             cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
00630             // Loop over all elements of the corresponding colvar value
00631             for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
00632                 // ds/dz, z = vector of CVs
00633                 tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
00634                 tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
00635                 // Apply the gradients to the atom groups in i-th cv
00636                 // Loop over all atom groups
00637                 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00638                     // Loop over all atoms in the k-th atom group
00639                     for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
00640                         // Chain rule
00641                         (*(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]);
00642                     }
00643                 }
00644             }
00645         }
00646     }
00647 }
00648 
00649 void colvar::gspathCV::apply_force(colvarvalue const &force) {
00650     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00651         // If this CV us explicit gradients, then atomic gradients is already calculated
00652         // We can apply the force to atom groups directly
00653         if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00654             for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00655                 (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
00656             }
00657         } else {
00658             // Temporary variables storing gradients
00659             colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
00660             colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
00661             // Compute factors for polynomial combinations
00662             cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
00663             for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
00664                 // ds/dz, z = vector of CVs
00665                 tmp_cv_grad_v1[j_elem] = -1.0 * sign * 0.5 * dfdv1[i_cv][j_elem] / M;
00666                 tmp_cv_grad_v2[j_elem] = sign * 0.5 * dfdv2[i_cv][j_elem] / M;
00667             }
00668             colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
00669             cv[i_cv]->apply_force(cv_force);
00670         }
00671     }
00672 }
00673 
00674 colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) {
00675     set_function_type("gzpathCV");
00676     cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n"));
00677     // Initialize variables for future calculation
00678     M = cvm::real(total_reference_frames - 1);
00679     m = 1.0;
00680     get_keyval(conf, "useSecondClosestFrame", use_second_closest_frame, true);
00681     if (use_second_closest_frame == true) {
00682         cvm::log(std::string("Geometric path z(σ) will use the second closest frame to compute s_(m-1)\n"));
00683     } else {
00684         cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m-1)\n"));
00685     }
00686     get_keyval(conf, "useThirdClosestFrame", use_third_closest_frame, false);
00687     if (use_third_closest_frame == true) {
00688         cvm::log(std::string("Geometric path z(σ) will use the third closest frame to compute s_(m+1)\n"));
00689     } else {
00690         cvm::log(std::string("Geometric path z(σ) will use the neighbouring frame to compute s_(m+1)\n"));
00691     }
00692     bool b_use_z_square = false;
00693     get_keyval(conf, "useZsquare", b_use_z_square, false);
00694     if (b_use_z_square == true) {
00695         cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n"));
00696     }
00697     if (total_reference_frames < 2) {
00698         cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpathCV requires at least 2 frames.\n");
00699         return;
00700     }
00701     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);
00702     x.type(colvarvalue::type_scalar);
00703 }
00704 
00705 colvar::gzpathCV::~gzpathCV() {
00706 }
00707 
00708 void colvar::gzpathCV::updateDistanceToReferenceFrames() {
00709     computeDistanceToReferenceFrames(frame_distances);
00710 }
00711 
00712 void colvar::gzpathCV::prepareVectors() {
00713     // Compute v1, v2 and v3
00714     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00715         // values of sub-cvc are computed in update_distances
00716         // cv[i_cv]->calc_value();
00717         colvarvalue f1_ref_cv_i_value(ref_cv[min_frame_index_1][i_cv]);
00718         colvarvalue f2_ref_cv_i_value(ref_cv[min_frame_index_2][i_cv]);
00719         colvarvalue current_cv_value(cv[i_cv]->value());
00720         // polynomial combination allowed
00721         if (current_cv_value.type() == colvarvalue::type_scalar) {
00722             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));
00723             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;
00724         } else {
00725             v1[i_cv] = f1_ref_cv_i_value - cv[i_cv]->sup_coeff * current_cv_value;
00726             v2[i_cv] = cv[i_cv]->sup_coeff * current_cv_value - f2_ref_cv_i_value;
00727         }
00728         v4[i_cv] = f1_ref_cv_i_value - f2_ref_cv_i_value;
00729         cv[i_cv]->wrap(v1[i_cv]);
00730         cv[i_cv]->wrap(v2[i_cv]);
00731         cv[i_cv]->wrap(v4[i_cv]);
00732     }
00733     if (min_frame_index_3 < 0 || min_frame_index_3 > M) {
00734         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00735             v3[i_cv] = ref_cv[min_frame_index_1][i_cv] - ref_cv[min_frame_index_2][i_cv];
00736             cv[i_cv]->wrap(v3[i_cv]);
00737         }
00738     } else {
00739         for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00740             v3[i_cv] = ref_cv[min_frame_index_3][i_cv] - ref_cv[min_frame_index_1][i_cv];
00741             cv[i_cv]->wrap(v3[i_cv]);
00742         }
00743     }
00744 }
00745 
00746 void colvar::gzpathCV::calc_value() {
00747     computeValue();
00748     x = z;
00749 }
00750 
00751 void colvar::gzpathCV::calc_gradients() {
00752     computeDerivatives();
00753     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00754         // No matter whether the i-th cv uses implicit gradient, compute it first.
00755         cv[i_cv]->calc_gradients();
00756         // If the gradient is not implicit, then add the gradients to its atom groups
00757         if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00758             // Temporary variables storing gradients
00759             colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
00760             colvarvalue tmp_cv_grad_v2 =  1.0 * dzdv2[i_cv];
00761             // Compute factors for polynomial combinations
00762             cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
00763             for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
00764                 // Apply the gradients to the atom groups in i-th cv
00765                 // Loop over all atom groups
00766                 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00767                     // Loop over all atoms in the k-th atom group
00768                     for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
00769                         // Chain rule
00770                         (*(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]);
00771                     }
00772                 }
00773             }
00774         }
00775     }
00776 }
00777 
00778 void colvar::gzpathCV::apply_force(colvarvalue const &force) {
00779     for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00780         // If this CV us explicit gradients, then atomic gradients is already calculated
00781         // We can apply the force to atom groups directly
00782         if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00783             for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00784                 (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value);
00785             }
00786         } else {
00787             colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
00788             colvarvalue tmp_cv_grad_v2 =  1.0 * dzdv2[i_cv];
00789             // Temporary variables storing gradients
00790             // Compute factors for polynomial combinations
00791             cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
00792             colvarvalue cv_force = force.real_value * factor_polynomial * (tmp_cv_grad_v1 + tmp_cv_grad_v2);
00793             cv[i_cv]->apply_force(cv_force);
00794         }
00795     }
00796 }
00797 
00798 #endif

Generated on Fri Nov 8 02:44:14 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002