00001 #if (__cplusplus >= 201103L)
00002
00003
00004
00005
00006
00007
00008
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
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
00032
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
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
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
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
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
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
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
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
00159 v1[i_atom] = reference_frames[min_frame_index_1][i_atom] - (*(comp_atoms[min_frame_index_1]))[i_atom].pos;
00160
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
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
00246
00247
00248
00249
00250
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
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
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
00342
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
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
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
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
00414 std::sort(cv.begin(), cv.end(), colvar::compare_cvc);
00415
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
00426
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
00494 if (current_cv_value.type() == colvarvalue::type_scalar) {
00495
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
00533
00534
00535
00536
00537
00538
00539
00540
00541 remove_all_children();
00542
00543
00544 for (auto it = cv.begin(); it != cv.end(); ++it) {
00545 delete (*it);
00546 }
00547
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
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
00583 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00584
00585
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
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
00622 cv[i_cv]->calc_gradients();
00623
00624 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00625
00626 colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
00627 colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
00628
00629 cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv);
00630
00631 for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) {
00632
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
00636
00637 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00638
00639 for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
00640
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
00652
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
00659 colvarvalue tmp_cv_grad_v1(cv[i_cv]->value());
00660 colvarvalue tmp_cv_grad_v2(cv[i_cv]->value());
00661
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
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
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
00714 for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) {
00715
00716
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
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
00755 cv[i_cv]->calc_gradients();
00756
00757 if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) {
00758
00759 colvarvalue tmp_cv_grad_v1 = -1.0 * dzdv1[i_cv];
00760 colvarvalue tmp_cv_grad_v2 = 1.0 * dzdv2[i_cv];
00761
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
00765
00766 for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) {
00767
00768 for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) {
00769
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
00781
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
00790
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