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

colvarcomp_distances.C

Go to the documentation of this file.
00001 // -*- c++ -*-
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 <algorithm>
00011 
00012 #include "colvarmodule.h"
00013 #include "colvarvalue.h"
00014 #include "colvarparse.h"
00015 #include "colvar.h"
00016 #include "colvarcomp.h"
00017 
00018 
00019 
00020 colvar::distance::distance(std::string const &conf)
00021   : cvc(conf)
00022 {
00023   set_function_type("distance");
00024   init_as_distance();
00025 
00026   provide(f_cvc_inv_gradient);
00027   provide(f_cvc_Jacobian);
00028   enable(f_cvc_com_based);
00029 
00030   group1 = parse_group(conf, "group1");
00031   group2 = parse_group(conf, "group2");
00032 
00033   init_total_force_params(conf);
00034 }
00035 
00036 
00037 colvar::distance::distance()
00038   : cvc()
00039 {
00040   set_function_type("distance");
00041   init_as_distance();
00042 
00043   provide(f_cvc_inv_gradient);
00044   provide(f_cvc_Jacobian);
00045   enable(f_cvc_com_based);
00046 }
00047 
00048 
00049 void colvar::distance::calc_value()
00050 {
00051   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00052     dist_v = group2->center_of_mass() - group1->center_of_mass();
00053   } else {
00054     dist_v = cvm::position_distance(group1->center_of_mass(),
00055                                     group2->center_of_mass());
00056   }
00057   x.real_value = dist_v.norm();
00058 }
00059 
00060 
00061 void colvar::distance::calc_gradients()
00062 {
00063   cvm::rvector const u = dist_v.unit();
00064   group1->set_weighted_gradient(-1.0 * u);
00065   group2->set_weighted_gradient(       u);
00066 }
00067 
00068 
00069 void colvar::distance::calc_force_invgrads()
00070 {
00071   group1->read_total_forces();
00072   if (is_enabled(f_cvc_one_site_total_force)) {
00073     ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
00074   } else {
00075     group2->read_total_forces();
00076     ft.real_value = 0.5 * ((group2->total_force() - group1->total_force()) * dist_v.unit());
00077   }
00078 }
00079 
00080 
00081 void colvar::distance::calc_Jacobian_derivative()
00082 {
00083   jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
00084 }
00085 
00086 
00087 void colvar::distance::apply_force(colvarvalue const &force)
00088 {
00089   if (!group1->noforce)
00090     group1->apply_colvar_force(force.real_value);
00091 
00092   if (!group2->noforce)
00093     group2->apply_colvar_force(force.real_value);
00094 }
00095 
00096 
00097 simple_scalar_dist_functions(distance)
00098 
00099 
00100 
00101 colvar::distance_vec::distance_vec(std::string const &conf)
00102   : distance(conf)
00103 {
00104   set_function_type("distanceVec");
00105   enable(f_cvc_com_based);
00106   disable(f_cvc_explicit_gradient);
00107   x.type(colvarvalue::type_3vector);
00108 }
00109 
00110 
00111 colvar::distance_vec::distance_vec()
00112   : distance()
00113 {
00114   set_function_type("distanceVec");
00115   enable(f_cvc_com_based);
00116   disable(f_cvc_explicit_gradient);
00117   x.type(colvarvalue::type_3vector);
00118 }
00119 
00120 
00121 void colvar::distance_vec::calc_value()
00122 {
00123   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00124     x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
00125   } else {
00126     x.rvector_value = cvm::position_distance(group1->center_of_mass(),
00127                                              group2->center_of_mass());
00128   }
00129 }
00130 
00131 
00132 void colvar::distance_vec::calc_gradients()
00133 {
00134   // gradients are not stored: a 3x3 matrix for each atom would be
00135   // needed to store just the identity matrix
00136 }
00137 
00138 
00139 void colvar::distance_vec::apply_force(colvarvalue const &force)
00140 {
00141   if (!group1->noforce)
00142     group1->apply_force(-1.0 * force.rvector_value);
00143 
00144   if (!group2->noforce)
00145     group2->apply_force(       force.rvector_value);
00146 }
00147 
00148 
00149 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
00150                                       colvarvalue const &x2) const
00151 {
00152   return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2();
00153 }
00154 
00155 
00156 colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
00157                                               colvarvalue const &x2) const
00158 {
00159   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
00160 }
00161 
00162 
00163 colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
00164                                               colvarvalue const &x2) const
00165 {
00166   return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
00167 }
00168 
00169 
00170 
00171 colvar::distance_z::distance_z(std::string const &conf)
00172   : cvc(conf)
00173 {
00174   set_function_type("distanceZ");
00175   provide(f_cvc_inv_gradient);
00176   provide(f_cvc_Jacobian);
00177   enable(f_cvc_com_based);
00178   x.type(colvarvalue::type_scalar);
00179 
00180   // TODO detect PBC from MD engine (in simple cases)
00181   // and then update period in real time
00182   if (period != 0.0) {
00183     enable(f_cvc_periodic);
00184   }
00185 
00186   if ((wrap_center != 0.0) && !is_enabled(f_cvc_periodic)) {
00187     cvm::error("Error: wrapAround was defined in a distanceZ component,"
00188                 " but its period has not been set.\n");
00189     return;
00190   }
00191 
00192   main = parse_group(conf, "main");
00193   ref1 = parse_group(conf, "ref");
00194   // this group is optional
00195   ref2 = parse_group(conf, "ref2", true);
00196 
00197   if ( ref2 ) {
00198     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"\n");
00199     fixed_axis = false;
00200     if (key_lookup(conf, "axis"))
00201       cvm::log("Warning: explicit axis definition will be ignored!\n");
00202   } else {
00203     if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) {
00204       if (axis.norm2() == 0.0) {
00205         cvm::error("Axis vector is zero!");
00206         return;
00207       }
00208       if (axis.norm2() != 1.0) {
00209         axis = axis.unit();
00210         cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n");
00211       }
00212     }
00213     fixed_axis = true;
00214   }
00215 
00216   init_total_force_params(conf);
00217 
00218 }
00219 
00220 
00221 colvar::distance_z::distance_z()
00222 {
00223   set_function_type("distanceZ");
00224   provide(f_cvc_inv_gradient);
00225   provide(f_cvc_Jacobian);
00226   enable(f_cvc_com_based);
00227   x.type(colvarvalue::type_scalar);
00228 }
00229 
00230 
00231 void colvar::distance_z::calc_value()
00232 {
00233   if (fixed_axis) {
00234     if (!is_enabled(f_cvc_pbc_minimum_image)) {
00235       dist_v = main->center_of_mass() - ref1->center_of_mass();
00236     } else {
00237       dist_v = cvm::position_distance(ref1->center_of_mass(),
00238                                       main->center_of_mass());
00239     }
00240   } else {
00241 
00242     if (!is_enabled(f_cvc_pbc_minimum_image)) {
00243       dist_v = main->center_of_mass() -
00244                (0.5 * (ref1->center_of_mass() + ref2->center_of_mass()));
00245       axis = ref2->center_of_mass() - ref1->center_of_mass();
00246     } else {
00247       dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() +
00248                                              ref2->center_of_mass()),
00249                                       main->center_of_mass());
00250       axis = cvm::position_distance(ref1->center_of_mass(),
00251                                     ref2->center_of_mass());
00252     }
00253     axis_norm = axis.norm();
00254     axis = axis.unit();
00255   }
00256   x.real_value = axis * dist_v;
00257   this->wrap(x);
00258 }
00259 
00260 
00261 void colvar::distance_z::calc_gradients()
00262 {
00263   main->set_weighted_gradient( axis );
00264 
00265   if (fixed_axis) {
00266     ref1->set_weighted_gradient(-1.0 * axis);
00267   } else {
00268     if (!is_enabled(f_cvc_pbc_minimum_image)) {
00269       ref1->set_weighted_gradient( 1.0 / axis_norm *
00270                                    (main->center_of_mass() - ref2->center_of_mass() -
00271                                    x.real_value * axis ));
00272       ref2->set_weighted_gradient( 1.0 / axis_norm *
00273                                    (ref1->center_of_mass() - main->center_of_mass() +
00274                                    x.real_value * axis ));
00275     } else {
00276       ref1->set_weighted_gradient( 1.0 / axis_norm * (
00277         cvm::position_distance(ref2->center_of_mass(),
00278                                main->center_of_mass()) - x.real_value * axis ));
00279       ref2->set_weighted_gradient( 1.0 / axis_norm * (
00280         cvm::position_distance(main->center_of_mass(),
00281                                ref1->center_of_mass()) + x.real_value * axis ));
00282     }
00283   }
00284 }
00285 
00286 
00287 void colvar::distance_z::calc_force_invgrads()
00288 {
00289   main->read_total_forces();
00290 
00291   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
00292     ref1->read_total_forces();
00293     ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
00294   } else {
00295     ft.real_value = main->total_force() * axis;
00296   }
00297 }
00298 
00299 
00300 void colvar::distance_z::calc_Jacobian_derivative()
00301 {
00302   jd.real_value = 0.0;
00303 }
00304 
00305 
00306 void colvar::distance_z::apply_force(colvarvalue const &force)
00307 {
00308   if (!ref1->noforce)
00309     ref1->apply_colvar_force(force.real_value);
00310 
00311   if (ref2 && !ref2->noforce)
00312     ref2->apply_colvar_force(force.real_value);
00313 
00314   if (!main->noforce)
00315     main->apply_colvar_force(force.real_value);
00316 }
00317 
00318 
00319 // Differences should always be wrapped around 0 (ignoring wrap_center)
00320 cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
00321                                     colvarvalue const &x2) const
00322 {
00323   cvm::real diff = x1.real_value - x2.real_value;
00324   if (is_enabled(f_cvc_periodic)) {
00325     cvm::real shift = cvm::floor(diff/period + 0.5);
00326     diff -= shift * period;
00327   }
00328   return diff * diff;
00329 }
00330 
00331 
00332 colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
00333                                             colvarvalue const &x2) const
00334 {
00335   cvm::real diff = x1.real_value - x2.real_value;
00336   if (is_enabled(f_cvc_periodic)) {
00337     cvm::real shift = cvm::floor(diff/period + 0.5);
00338     diff -= shift * period;
00339   }
00340   return 2.0 * diff;
00341 }
00342 
00343 
00344 colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
00345                                             colvarvalue const &x2) const
00346 {
00347   cvm::real diff = x1.real_value - x2.real_value;
00348   if (is_enabled(f_cvc_periodic)) {
00349     cvm::real shift = cvm::floor(diff/period + 0.5);
00350     diff -= shift * period;
00351   }
00352   return (-2.0) * diff;
00353 }
00354 
00355 
00356 void colvar::distance_z::wrap(colvarvalue &x_unwrapped) const
00357 {
00358   if (!is_enabled(f_cvc_periodic)) {
00359     // don't wrap if the period has not been set
00360     return;
00361   }
00362   cvm::real shift =
00363     cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5);
00364   x_unwrapped.real_value -= shift * period;
00365 }
00366 
00367 
00368 
00369 colvar::distance_xy::distance_xy(std::string const &conf)
00370   : distance_z(conf)
00371 {
00372   set_function_type("distanceXY");
00373   init_as_distance();
00374 
00375   provide(f_cvc_inv_gradient);
00376   provide(f_cvc_Jacobian);
00377   enable(f_cvc_com_based);
00378 }
00379 
00380 
00381 colvar::distance_xy::distance_xy()
00382   : distance_z()
00383 {
00384   set_function_type("distanceXY");
00385   init_as_distance();
00386 
00387   provide(f_cvc_inv_gradient);
00388   provide(f_cvc_Jacobian);
00389   enable(f_cvc_com_based);
00390 }
00391 
00392 
00393 void colvar::distance_xy::calc_value()
00394 {
00395   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00396     dist_v = main->center_of_mass() - ref1->center_of_mass();
00397   } else {
00398     dist_v = cvm::position_distance(ref1->center_of_mass(),
00399                                     main->center_of_mass());
00400   }
00401   if (!fixed_axis) {
00402     if (!is_enabled(f_cvc_pbc_minimum_image)) {
00403       v12 = ref2->center_of_mass() - ref1->center_of_mass();
00404     } else {
00405       v12 = cvm::position_distance(ref1->center_of_mass(),
00406                                    ref2->center_of_mass());
00407     }
00408     axis_norm = v12.norm();
00409     axis = v12.unit();
00410   }
00411 
00412   dist_v_ortho = dist_v - (dist_v * axis) * axis;
00413   x.real_value = dist_v_ortho.norm();
00414 }
00415 
00416 
00417 void colvar::distance_xy::calc_gradients()
00418 {
00419   // Intermediate quantity (r_P3 / r_12 where P is the projection
00420   // of 3(main) on the plane orthogonal to 12, containing 1 (ref1))
00421   cvm::real A;
00422   cvm::real x_inv;
00423 
00424   if (x.real_value == 0.0) return;
00425   x_inv = 1.0 / x.real_value;
00426 
00427   if (fixed_axis) {
00428     ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);
00429     main->set_weighted_gradient(       x_inv * dist_v_ortho);
00430   } else {
00431     if (!is_enabled(f_cvc_pbc_minimum_image)) {
00432       v13 = main->center_of_mass() - ref1->center_of_mass();
00433     } else {
00434       v13 = cvm::position_distance(ref1->center_of_mass(),
00435                                    main->center_of_mass());
00436     }
00437     A = (dist_v * axis) / axis_norm;
00438 
00439     ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho);
00440     ref2->set_weighted_gradient( -A        * x_inv * dist_v_ortho);
00441     main->set_weighted_gradient(      1.0  * x_inv * dist_v_ortho);
00442   }
00443 }
00444 
00445 
00446 void colvar::distance_xy::calc_force_invgrads()
00447 {
00448   main->read_total_forces();
00449 
00450   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
00451     ref1->read_total_forces();
00452     ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
00453   } else {
00454     ft.real_value = 1.0 / x.real_value * main->total_force() * dist_v_ortho;
00455   }
00456 }
00457 
00458 
00459 void colvar::distance_xy::calc_Jacobian_derivative()
00460 {
00461   jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
00462 }
00463 
00464 
00465 void colvar::distance_xy::apply_force(colvarvalue const &force)
00466 {
00467   if (!ref1->noforce)
00468     ref1->apply_colvar_force(force.real_value);
00469 
00470   if (ref2 && !ref2->noforce)
00471     ref2->apply_colvar_force(force.real_value);
00472 
00473   if (!main->noforce)
00474     main->apply_colvar_force(force.real_value);
00475 }
00476 
00477 
00478 simple_scalar_dist_functions(distance_xy)
00479 
00480 
00481 
00482 colvar::distance_dir::distance_dir(std::string const &conf)
00483   : distance(conf)
00484 {
00485   set_function_type("distanceDir");
00486   enable(f_cvc_com_based);
00487   disable(f_cvc_explicit_gradient);
00488   x.type(colvarvalue::type_unit3vector);
00489 }
00490 
00491 
00492 colvar::distance_dir::distance_dir()
00493   : distance()
00494 {
00495   set_function_type("distanceDir");
00496   enable(f_cvc_com_based);
00497   disable(f_cvc_explicit_gradient);
00498   x.type(colvarvalue::type_unit3vector);
00499 }
00500 
00501 
00502 void colvar::distance_dir::calc_value()
00503 {
00504   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00505     dist_v = group2->center_of_mass() - group1->center_of_mass();
00506   } else {
00507     dist_v = cvm::position_distance(group1->center_of_mass(),
00508                                     group2->center_of_mass());
00509   }
00510   x.rvector_value = dist_v.unit();
00511 }
00512 
00513 
00514 void colvar::distance_dir::calc_gradients()
00515 {
00516   // gradients are computed on the fly within apply_force()
00517   // Note: could be a problem if a future bias relies on gradient
00518   // calculations...
00519   // TODO in new deps system: remove dependency of biasing force to gradient?
00520   // That way we could tell apart an explicit gradient dependency
00521 }
00522 
00523 
00524 void colvar::distance_dir::apply_force(colvarvalue const &force)
00525 {
00526   // remove the radial force component
00527   cvm::real const iprod = force.rvector_value * x.rvector_value;
00528   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
00529 
00530   if (!group1->noforce)
00531     group1->apply_force(-1.0 * force_tang);
00532 
00533   if (!group2->noforce)
00534     group2->apply_force(       force_tang);
00535 }
00536 
00537 
00538 cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
00539                                       colvarvalue const &x2) const
00540 {
00541   return (x1.rvector_value - x2.rvector_value).norm2();
00542 }
00543 
00544 
00545 colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
00546                                               colvarvalue const &x2) const
00547 {
00548   return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vectorderiv);
00549 }
00550 
00551 
00552 colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
00553                                               colvarvalue const &x2) const
00554 {
00555   return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vectorderiv);
00556 }
00557 
00558 
00559 
00560 colvar::distance_inv::distance_inv(std::string const &conf)
00561   : cvc(conf)
00562 {
00563   set_function_type("distanceInv");
00564   init_as_distance();
00565 
00566   group1 = parse_group(conf, "group1");
00567   group2 = parse_group(conf, "group2");
00568 
00569   get_keyval(conf, "exponent", exponent, 6);
00570   if (exponent%2) {
00571     cvm::error("Error: odd exponent provided, can only use even ones.\n");
00572     return;
00573   }
00574   if (exponent <= 0) {
00575     cvm::error("Error: negative or zero exponent provided.\n");
00576     return;
00577   }
00578 
00579   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00580     for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
00581       if (ai1->id == ai2->id) {
00582         cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n");
00583         return;
00584       }
00585     }
00586   }
00587 
00588   if (is_enabled(f_cvc_debug_gradient)) {
00589     cvm::log("Warning: debugGradients will not give correct results "
00590              "for distanceInv, because its value and gradients are computed "
00591              "simultaneously.\n");
00592   }
00593 }
00594 
00595 
00596 void colvar::distance_inv::calc_value()
00597 {
00598   x.real_value = 0.0;
00599   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00600     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00601       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
00602         cvm::rvector const dv = ai2->pos - ai1->pos;
00603         cvm::real const d2 = dv.norm2();
00604         cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
00605         x.real_value += dinv;
00606         cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
00607         ai1->grad += -1.0 * dsumddv;
00608         ai2->grad +=        dsumddv;
00609       }
00610     }
00611   } else {
00612     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00613       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
00614         cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos);
00615         cvm::real const d2 = dv.norm2();
00616         cvm::real const dinv = cvm::integer_power(d2, -1*(exponent/2));
00617         x.real_value += dinv;
00618         cvm::rvector const dsumddv = -1.0*(exponent/2) * dinv/d2 * 2.0 * dv;
00619         ai1->grad += -1.0 * dsumddv;
00620         ai2->grad +=        dsumddv;
00621       }
00622     }
00623   }
00624 
00625   x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
00626   x.real_value = cvm::pow(x.real_value, -1.0/cvm::real(exponent));
00627 
00628   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) *
00629     cvm::integer_power(x.real_value, exponent+1) /
00630     cvm::real(group1->size() * group2->size());
00631   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00632     ai1->grad *= dxdsum;
00633   }
00634   for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
00635     ai2->grad *= dxdsum;
00636   }
00637 }
00638 
00639 
00640 void colvar::distance_inv::calc_gradients()
00641 {
00642 }
00643 
00644 
00645 void colvar::distance_inv::apply_force(colvarvalue const &force)
00646 {
00647   if (!group1->noforce)
00648     group1->apply_colvar_force(force.real_value);
00649 
00650   if (!group2->noforce)
00651     group2->apply_colvar_force(force.real_value);
00652 }
00653 
00654 
00655 simple_scalar_dist_functions(distance_inv)
00656 
00657 
00658 
00659 colvar::distance_pairs::distance_pairs(std::string const &conf)
00660   : cvc(conf)
00661 {
00662   set_function_type("distancePairs");
00663 
00664   group1 = parse_group(conf, "group1");
00665   group2 = parse_group(conf, "group2");
00666 
00667   x.type(colvarvalue::type_vector);
00668   disable(f_cvc_explicit_gradient);
00669   x.vector1d_value.resize(group1->size() * group2->size());
00670 }
00671 
00672 
00673 colvar::distance_pairs::distance_pairs()
00674 {
00675   set_function_type("distancePairs");
00676   disable(f_cvc_explicit_gradient);
00677   x.type(colvarvalue::type_vector);
00678 }
00679 
00680 
00681 void colvar::distance_pairs::calc_value()
00682 {
00683   x.vector1d_value.resize(group1->size() * group2->size());
00684 
00685   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00686     size_t i1, i2;
00687     for (i1 = 0; i1 < group1->size(); i1++) {
00688       for (i2 = 0; i2 < group2->size(); i2++) {
00689         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
00690         cvm::real const d = dv.norm();
00691         x.vector1d_value[i1*group2->size() + i2] = d;
00692         (*group1)[i1].grad = -1.0 * dv.unit();
00693         (*group2)[i2].grad =  dv.unit();
00694       }
00695     }
00696   } else {
00697     size_t i1, i2;
00698     for (i1 = 0; i1 < group1->size(); i1++) {
00699       for (i2 = 0; i2 < group2->size(); i2++) {
00700         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
00701                                                        (*group2)[i2].pos);
00702         cvm::real const d = dv.norm();
00703         x.vector1d_value[i1*group2->size() + i2] = d;
00704         (*group1)[i1].grad = -1.0 * dv.unit();
00705         (*group2)[i2].grad =  dv.unit();
00706       }
00707     }
00708   }
00709 }
00710 
00711 
00712 void colvar::distance_pairs::calc_gradients()
00713 {
00714   // will be calculated on the fly in apply_force()
00715 }
00716 
00717 
00718 void colvar::distance_pairs::apply_force(colvarvalue const &force)
00719 {
00720   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00721     size_t i1, i2;
00722     for (i1 = 0; i1 < group1->size(); i1++) {
00723       for (i2 = 0; i2 < group2->size(); i2++) {
00724         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
00725         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
00726         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
00727       }
00728     }
00729   } else {
00730     size_t i1, i2;
00731     for (i1 = 0; i1 < group1->size(); i1++) {
00732       for (i2 = 0; i2 < group2->size(); i2++) {
00733         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos,
00734                                                        (*group2)[i2].pos);
00735         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
00736         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
00737       }
00738     }
00739   }
00740 }
00741 
00742 
00743 
00744 colvar::dipole_magnitude::dipole_magnitude(std::string const &conf)
00745   : cvc(conf)
00746 {
00747   set_function_type("dipoleMagnitude");
00748   atoms = parse_group(conf, "atoms");
00749   init_total_force_params(conf);
00750   x.type(colvarvalue::type_scalar);
00751 }
00752 
00753 
00754 colvar::dipole_magnitude::dipole_magnitude(cvm::atom const &a1)
00755 {
00756   set_function_type("dipoleMagnitude");
00757   atoms = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
00758   register_atom_group(atoms);
00759   x.type(colvarvalue::type_scalar);
00760 }
00761 
00762 
00763 colvar::dipole_magnitude::dipole_magnitude()
00764 {
00765   set_function_type("dipoleMagnitude");
00766   x.type(colvarvalue::type_scalar);
00767 }
00768 
00769 
00770 void colvar::dipole_magnitude::calc_value()
00771 {
00772   cvm::atom_pos const atomsCom = atoms->center_of_mass();
00773   atoms->calc_dipole(atomsCom);
00774   dipoleV = atoms->dipole();
00775   x.real_value = dipoleV.norm();
00776 }
00777 
00778 
00779 void colvar::dipole_magnitude::calc_gradients()
00780 {
00781   cvm::real const aux1 = atoms->total_charge/atoms->total_mass;
00782   cvm::atom_pos const dipVunit = dipoleV.unit();
00783 
00784   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00785     ai->grad = (ai->charge - aux1*ai->mass) * dipVunit;
00786   }
00787 }
00788 
00789 
00790 void colvar::dipole_magnitude::apply_force(colvarvalue const &force)
00791 {
00792   if (!atoms->noforce) {
00793     atoms->apply_colvar_force(force.real_value);
00794   }
00795 }
00796 
00797 
00798 simple_scalar_dist_functions(dipole_magnitude)
00799 
00800 
00801 
00802 colvar::gyration::gyration(std::string const &conf)
00803   : cvc(conf)
00804 {
00805   set_function_type("gyration");
00806   init_as_distance();
00807 
00808   provide(f_cvc_inv_gradient);
00809   provide(f_cvc_Jacobian);
00810   atoms = parse_group(conf, "atoms");
00811 
00812   if (atoms->b_user_defined_fit) {
00813     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
00814   } else {
00815     atoms->enable(f_ag_center);
00816     atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));
00817     atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0));
00818   }
00819 }
00820 
00821 
00822 void colvar::gyration::calc_value()
00823 {
00824   x.real_value = 0.0;
00825   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00826     x.real_value += (ai->pos).norm2();
00827   }
00828   x.real_value = cvm::sqrt(x.real_value / cvm::real(atoms->size()));
00829 }
00830 
00831 
00832 void colvar::gyration::calc_gradients()
00833 {
00834   cvm::real const drdx = 1.0/(cvm::real(atoms->size()) * x.real_value);
00835   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00836     ai->grad = drdx * ai->pos;
00837   }
00838 }
00839 
00840 
00841 void colvar::gyration::calc_force_invgrads()
00842 {
00843   atoms->read_total_forces();
00844 
00845   cvm::real const dxdr = 1.0/x.real_value;
00846   ft.real_value = 0.0;
00847 
00848   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00849     ft.real_value += dxdr * ai->pos * ai->total_force;
00850   }
00851 }
00852 
00853 
00854 void colvar::gyration::calc_Jacobian_derivative()
00855 {
00856   jd = x.real_value ? (3.0 * cvm::real(atoms->size()) - 4.0) / x.real_value : 0.0;
00857 }
00858 
00859 
00860 void colvar::gyration::apply_force(colvarvalue const &force)
00861 {
00862   if (!atoms->noforce)
00863     atoms->apply_colvar_force(force.real_value);
00864 }
00865 
00866 
00867 simple_scalar_dist_functions(gyration)
00868 
00869 
00870 
00871 colvar::inertia::inertia(std::string const &conf)
00872   : gyration(conf)
00873 {
00874   set_function_type("inertia");
00875   init_as_distance();
00876 }
00877 
00878 
00879 void colvar::inertia::calc_value()
00880 {
00881   x.real_value = 0.0;
00882   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00883     x.real_value += (ai->pos).norm2();
00884   }
00885 }
00886 
00887 
00888 void colvar::inertia::calc_gradients()
00889 {
00890   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00891     ai->grad = 2.0 * ai->pos;
00892   }
00893 }
00894 
00895 
00896 void colvar::inertia::apply_force(colvarvalue const &force)
00897 {
00898   if (!atoms->noforce)
00899     atoms->apply_colvar_force(force.real_value);
00900 }
00901 
00902 
00903 simple_scalar_dist_functions(inertia_z)
00904 
00905 
00906 
00907 colvar::inertia_z::inertia_z(std::string const &conf)
00908   : inertia(conf)
00909 {
00910   set_function_type("inertiaZ");
00911   init_as_distance();
00912   if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) {
00913     if (axis.norm2() == 0.0) {
00914       cvm::error("Axis vector is zero!", COLVARS_INPUT_ERROR);
00915       return;
00916     }
00917     if (axis.norm2() != 1.0) {
00918       axis = axis.unit();
00919       cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n");
00920     }
00921   }
00922 }
00923 
00924 
00925 void colvar::inertia_z::calc_value()
00926 {
00927   x.real_value = 0.0;
00928   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00929     cvm::real const iprod = ai->pos * axis;
00930     x.real_value += iprod * iprod;
00931   }
00932 }
00933 
00934 
00935 void colvar::inertia_z::calc_gradients()
00936 {
00937   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
00938     ai->grad = 2.0 * (ai->pos * axis) * axis;
00939   }
00940 }
00941 
00942 
00943 void colvar::inertia_z::apply_force(colvarvalue const &force)
00944 {
00945   if (!atoms->noforce)
00946     atoms->apply_colvar_force(force.real_value);
00947 }
00948 
00949 
00950 simple_scalar_dist_functions(inertia)
00951 
00952 
00953 
00954 
00955 colvar::rmsd::rmsd(std::string const &conf)
00956   : cvc(conf)
00957 {
00958   set_function_type("rmsd");
00959   init_as_distance();
00960 
00961   provide(f_cvc_inv_gradient);
00962 
00963   atoms = parse_group(conf, "atoms");
00964 
00965   if (!atoms || atoms->size() == 0) {
00966     cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD.");
00967     return;
00968   }
00969 
00970   bool b_Jacobian_derivative = true;
00971   if (atoms->fitting_group != NULL && b_Jacobian_derivative) {
00972     cvm::log("The option \"fittingGroup\" (alternative group for fitting) was enabled: "
00973               "Jacobian derivatives of the RMSD will not be calculated.\n");
00974     b_Jacobian_derivative = false;
00975   }
00976   if (b_Jacobian_derivative) provide(f_cvc_Jacobian);
00977 
00978   // the following is a simplified version of the corresponding atom group options;
00979   // we need this because the reference coordinates defined inside the atom group
00980   // may be used only for fitting, and even more so if fitting_group is used
00981   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
00982     cvm::log("Using reference positions from configuration file to calculate the variable.\n");
00983     if (ref_pos.size() != atoms->size()) {
00984       cvm::error("Error: the number of reference positions provided ("+
00985                   cvm::to_str(ref_pos.size())+
00986                   ") does not match the number of atoms of group \"atoms\" ("+
00987                   cvm::to_str(atoms->size())+").\n");
00988       return;
00989     }
00990   } else { // Only look for ref pos file if ref positions not already provided
00991     std::string ref_pos_file;
00992     if (get_keyval(conf, "refPositionsFile", ref_pos_file, std::string(""))) {
00993 
00994       if (ref_pos.size()) {
00995         cvm::error("Error: cannot specify \"refPositionsFile\" and "
00996                           "\"refPositions\" at the same time.\n");
00997         return;
00998       }
00999 
01000       std::string ref_pos_col;
01001       double ref_pos_col_value=0.0;
01002 
01003       if (get_keyval(conf, "refPositionsCol", ref_pos_col, std::string(""))) {
01004         // if provided, use PDB column to select coordinates
01005         bool found = get_keyval(conf, "refPositionsColValue", ref_pos_col_value, 0.0);
01006         if (found && ref_pos_col_value==0.0) {
01007           cvm::error("Error: refPositionsColValue, "
01008                      "if provided, must be non-zero.\n");
01009           return;
01010         }
01011       }
01012 
01013       ref_pos.resize(atoms->size());
01014 
01015       cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms,
01016                        ref_pos_col, ref_pos_col_value);
01017     } else {
01018       cvm::error("Error: no reference positions for RMSD; use either refPositions of refPositionsFile.");
01019       return;
01020     }
01021   }
01022 
01023   if (ref_pos.size() != atoms->size()) {
01024     cvm::error("Error: found " + cvm::to_str(ref_pos.size()) +
01025                     " reference positions for RMSD; expected " + cvm::to_str(atoms->size()));
01026     return;
01027   }
01028 
01029   if (atoms->b_user_defined_fit) {
01030     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
01031   } else {
01032     // Default: fit everything
01033     cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating it as a variable: "
01034               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
01035     atoms->enable(f_ag_center);
01036     atoms->enable(f_ag_rotate);
01037     // default case: reference positions for calculating the rmsd are also those used
01038     // for fitting
01039     atoms->ref_pos = ref_pos;
01040     atoms->center_ref_pos();
01041 
01042     cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
01043               "will not be computed as they cancel out in the gradients.");
01044     atoms->disable(f_ag_fit_gradients);
01045 
01046     // request the calculation of the derivatives of the rotation defined by the atom group
01047     atoms->rot.request_group1_gradients(atoms->size());
01048     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
01049     // this is only required for ABF, but we do both groups here for better caching
01050     atoms->rot.request_group2_gradients(atoms->size());
01051   }
01052 
01053   std::string perm_conf;
01054   size_t pos = 0; // current position in config string
01055   n_permutations = 1;
01056 
01057   while (key_lookup(conf, "atomPermutation", &perm_conf, &pos)) {
01058     cvm::main()->cite_feature("Symmetry-adapted RMSD");
01059     std::vector<size_t> perm;
01060     if (perm_conf.size()) {
01061       std::istringstream is(perm_conf);
01062       size_t index;
01063       while (is >> index) {
01064         std::vector<int> const &ids = atoms->ids();
01065         size_t const ia = std::find(ids.begin(), ids.end(), index-1) - ids.begin();
01066         if (ia == atoms->size()) {
01067           cvm::error("Error: atom id " + cvm::to_str(index) +
01068                     " is not a member of group \"atoms\".");
01069           return;
01070         }
01071         if (std::find(perm.begin(), perm.end(), ia) != perm.end()) {
01072           cvm::error("Error: atom id " + cvm::to_str(index) +
01073                     " is mentioned more than once in atomPermutation list.");
01074           return;
01075         }
01076         perm.push_back(ia);
01077       }
01078       if (perm.size() != atoms->size()) {
01079         cvm::error("Error: symmetry permutation in input contains " + cvm::to_str(perm.size()) +
01080                   " indices, but group \"atoms\" contains " + cvm::to_str(atoms->size()) + " atoms.");
01081         return;
01082       }
01083       cvm::log("atomPermutation = " + cvm::to_str(perm));
01084       n_permutations++;
01085       // Record a copy of reference positions in new order
01086       for (size_t ia = 0; ia < atoms->size(); ia++) {
01087         ref_pos.push_back(ref_pos[perm[ia]]);
01088       }
01089     }
01090   }
01091 }
01092 
01093 
01094 void colvar::rmsd::calc_value()
01095 {
01096   // rotational-translational fit is handled by the atom group
01097 
01098   x.real_value = 0.0;
01099   for (size_t ia = 0; ia < atoms->size(); ia++) {
01100     x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2();
01101   }
01102   best_perm_index = 0;
01103 
01104   // Compute sum of squares for each symmetry permutation of atoms, keep the smallest
01105   size_t ref_pos_index = atoms->size();
01106   for (size_t ip = 1; ip < n_permutations; ip++) {
01107     cvm::real value = 0.0;
01108     for (size_t ia = 0; ia < atoms->size(); ia++) {
01109       value += ((*atoms)[ia].pos - ref_pos[ref_pos_index++]).norm2();
01110     }
01111     if (value < x.real_value) {
01112       x.real_value = value;
01113       best_perm_index = ip;
01114     }
01115   }
01116   x.real_value /= cvm::real(atoms->size()); // MSD
01117   x.real_value = cvm::sqrt(x.real_value);
01118 }
01119 
01120 
01121 void colvar::rmsd::calc_gradients()
01122 {
01123   cvm::real const drmsddx2 = (x.real_value > 0.0) ?
01124     0.5 / (x.real_value * cvm::real(atoms->size())) :
01125     0.0;
01126 
01127   // Use the appropriate symmetry permutation of reference positions to calculate gradients
01128   size_t const start = atoms->size() * best_perm_index;
01129   for (size_t ia = 0; ia < atoms->size(); ia++) {
01130     (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[start + ia]));
01131   }
01132 }
01133 
01134 
01135 void colvar::rmsd::apply_force(colvarvalue const &force)
01136 {
01137   if (!atoms->noforce)
01138     atoms->apply_colvar_force(force.real_value);
01139 }
01140 
01141 
01142 void colvar::rmsd::calc_force_invgrads()
01143 {
01144   atoms->read_total_forces();
01145   ft.real_value = 0.0;
01146 
01147   // Note: gradient square norm is 1/N_atoms
01148 
01149   for (size_t ia = 0; ia < atoms->size(); ia++) {
01150     ft.real_value += (*atoms)[ia].grad * (*atoms)[ia].total_force;
01151   }
01152   ft.real_value *= atoms->size();
01153 }
01154 
01155 
01156 void colvar::rmsd::calc_Jacobian_derivative()
01157 {
01158   // divergence of the rotated coordinates (including only derivatives of the rotation matrix)
01159   cvm::real rotation_term = 0.0;
01160 
01161   // The rotation term only applies is coordinates are rotated
01162   if (atoms->is_enabled(f_ag_rotate)) {
01163 
01164     // gradient of the rotation matrix
01165     cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
01166     // gradients of products of 2 quaternion components
01167     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
01168     for (size_t ia = 0; ia < atoms->size(); ia++) {
01169 
01170       // Gradient of optimal quaternion wrt current Cartesian position
01171       cvm::vector1d<cvm::rvector> &dq = atoms->rot.dQ0_1[ia];
01172 
01173       g11 = 2.0 * (atoms->rot.q)[1]*dq[1];
01174       g22 = 2.0 * (atoms->rot.q)[2]*dq[2];
01175       g33 = 2.0 * (atoms->rot.q)[3]*dq[3];
01176       g01 = (atoms->rot.q)[0]*dq[1] + (atoms->rot.q)[1]*dq[0];
01177       g02 = (atoms->rot.q)[0]*dq[2] + (atoms->rot.q)[2]*dq[0];
01178       g03 = (atoms->rot.q)[0]*dq[3] + (atoms->rot.q)[3]*dq[0];
01179       g12 = (atoms->rot.q)[1]*dq[2] + (atoms->rot.q)[2]*dq[1];
01180       g13 = (atoms->rot.q)[1]*dq[3] + (atoms->rot.q)[3]*dq[1];
01181       g23 = (atoms->rot.q)[2]*dq[3] + (atoms->rot.q)[3]*dq[2];
01182 
01183       // Gradient of the rotation matrix wrt current Cartesian position
01184       grad_rot_mat[0][0] = -2.0 * (g22 + g33);
01185       grad_rot_mat[1][0] =  2.0 * (g12 + g03);
01186       grad_rot_mat[2][0] =  2.0 * (g13 - g02);
01187       grad_rot_mat[0][1] =  2.0 * (g12 - g03);
01188       grad_rot_mat[1][1] = -2.0 * (g11 + g33);
01189       grad_rot_mat[2][1] =  2.0 * (g01 + g23);
01190       grad_rot_mat[0][2] =  2.0 * (g02 + g13);
01191       grad_rot_mat[1][2] =  2.0 * (g23 - g01);
01192       grad_rot_mat[2][2] = -2.0 * (g11 + g22);
01193 
01194       cvm::atom_pos &y = ref_pos[ia];
01195 
01196       for (size_t alpha = 0; alpha < 3; alpha++) {
01197         for (size_t beta = 0; beta < 3; beta++) {
01198           rotation_term += grad_rot_mat[beta][alpha][alpha] * y[beta];
01199         // Note: equation was derived for inverse rotation (see colvars paper)
01200         // so here the matrix is transposed
01201         // (eq would give   divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];)
01202         }
01203       }
01204     }
01205   }
01206 
01207   // The translation term only applies is coordinates are centered
01208   cvm::real translation_term = atoms->is_enabled(f_ag_center) ? 3.0 : 0.0;
01209 
01210   jd.real_value = x.real_value > 0.0 ?
01211     (3.0 * atoms->size() - 1.0 - translation_term - rotation_term) / x.real_value :
01212     0.0;
01213 }
01214 
01215 
01216 simple_scalar_dist_functions(rmsd)
01217 
01218 
01219 
01220 colvar::eigenvector::eigenvector(std::string const &conf)
01221   : cvc(conf)
01222 {
01223   set_function_type("eigenvector");
01224   provide(f_cvc_inv_gradient);
01225   provide(f_cvc_Jacobian);
01226   x.type(colvarvalue::type_scalar);
01227 
01228   atoms = parse_group(conf, "atoms");
01229 
01230   {
01231     bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos);
01232 
01233     if (b_inline) {
01234       cvm::log("Using reference positions from input file.\n");
01235       if (ref_pos.size() != atoms->size()) {
01236         cvm::error("Error: reference positions do not "
01237                    "match the number of requested atoms.\n");
01238         return;
01239       }
01240     }
01241 
01242     std::string file_name;
01243     if (get_keyval(conf, "refPositionsFile", file_name)) {
01244 
01245       if (b_inline) {
01246         cvm::error("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
01247         return;
01248       }
01249 
01250       std::string file_col;
01251       double file_col_value=0.0;
01252       if (get_keyval(conf, "refPositionsCol", file_col, std::string(""))) {
01253         // use PDB flags if column is provided
01254         bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0);
01255         if (found && file_col_value==0.0) {
01256           cvm::error("Error: refPositionsColValue, "
01257                             "if provided, must be non-zero.\n");
01258           return;
01259         }
01260       }
01261 
01262       ref_pos.resize(atoms->size());
01263       cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
01264                        file_col, file_col_value);
01265     }
01266   }
01267 
01268   if (ref_pos.size() == 0) {
01269     cvm::error("Error: reference positions were not provided.\n", COLVARS_INPUT_ERROR);
01270     return;
01271   }
01272 
01273   if (ref_pos.size() != atoms->size()) {
01274     cvm::error("Error: reference positions do not "
01275                "match the number of requested atoms.\n", COLVARS_INPUT_ERROR);
01276     return;
01277   }
01278 
01279   // save for later the geometric center of the provided positions (may not be the origin)
01280   cvm::rvector ref_pos_center(0.0, 0.0, 0.0);
01281   for (size_t i = 0; i < atoms->size(); i++) {
01282     ref_pos_center += ref_pos[i];
01283   }
01284   ref_pos_center *= 1.0 / atoms->size();
01285 
01286   if (atoms->b_user_defined_fit) {
01287     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
01288   } else {
01289     // default: fit everything
01290     cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating the vector projection: "
01291               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
01292     atoms->enable(f_ag_center);
01293     atoms->enable(f_ag_rotate);
01294     atoms->ref_pos = ref_pos;
01295     atoms->center_ref_pos();
01296     atoms->disable(f_ag_fit_gradients); // cancel out if group is fitted on itself
01297                                         // and cvc is translationally invariant
01298 
01299     // request the calculation of the derivatives of the rotation defined by the atom group
01300     atoms->rot.request_group1_gradients(atoms->size());
01301     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
01302     // this is only required for ABF, but we do both groups here for better caching
01303     atoms->rot.request_group2_gradients(atoms->size());
01304   }
01305 
01306   {
01307     bool const b_inline = get_keyval(conf, "vector", eigenvec, eigenvec);
01308     // now load the eigenvector
01309     if (b_inline) {
01310       cvm::log("Using vector components from input file.\n");
01311       if (eigenvec.size() != atoms->size()) {
01312         cvm::error("Error: vector components do not "
01313                           "match the number of requested atoms->\n");
01314         return;
01315       }
01316     }
01317 
01318     std::string file_name;
01319     if (get_keyval(conf, "vectorFile", file_name)) {
01320 
01321       if (b_inline) {
01322         cvm::error("Error: vector and vectorFile cannot be specified at the same time.\n");
01323         return;
01324       }
01325 
01326       std::string file_col;
01327       double file_col_value=0.0;
01328       if (get_keyval(conf, "vectorCol", file_col, std::string(""))) {
01329         // use PDB flags if column is provided
01330         bool found = get_keyval(conf, "vectorColValue", file_col_value, 0.0);
01331         if (found && file_col_value==0.0) {
01332           cvm::error("Error: vectorColValue, if provided, must be non-zero.\n");
01333           return;
01334         }
01335       }
01336 
01337       eigenvec.resize(atoms->size());
01338       cvm::load_coords(file_name.c_str(), &eigenvec, atoms,
01339                        file_col, file_col_value);
01340     }
01341   }
01342 
01343   if (!ref_pos.size() || !eigenvec.size()) {
01344     cvm::error("Error: both reference coordinates"
01345                       "and eigenvector must be defined.\n");
01346     return;
01347   }
01348 
01349   cvm::atom_pos eig_center(0.0, 0.0, 0.0);
01350   for (size_t eil = 0; eil < atoms->size(); eil++) {
01351     eig_center += eigenvec[eil];
01352   }
01353   eig_center *= 1.0 / atoms->size();
01354   cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n");
01355 
01356   bool b_difference_vector = false;
01357   get_keyval(conf, "differenceVector", b_difference_vector, false);
01358 
01359   if (b_difference_vector) {
01360 
01361     if (atoms->is_enabled(f_ag_center)) {
01362       // both sets should be centered on the origin for fitting
01363       for (size_t i = 0; i < atoms->size(); i++) {
01364         eigenvec[i] -= eig_center;
01365         ref_pos[i]  -= ref_pos_center;
01366       }
01367     }
01368     if (atoms->is_enabled(f_ag_rotate)) {
01369       atoms->rot.calc_optimal_rotation(eigenvec, ref_pos);
01370       for (size_t i = 0; i < atoms->size(); i++) {
01371         eigenvec[i] = atoms->rot.rotate(eigenvec[i]);
01372       }
01373     }
01374     cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = x_vec - x_ref.\n");
01375     for (size_t i = 0; i < atoms->size(); i++) {
01376       eigenvec[i] -= ref_pos[i];
01377     }
01378     if (atoms->is_enabled(f_ag_center)) {
01379       // bring back the ref positions to where they were
01380       for (size_t i = 0; i < atoms->size(); i++) {
01381         ref_pos[i] += ref_pos_center;
01382       }
01383     }
01384 
01385   } else {
01386     cvm::log("Centering the provided vector to zero.\n");
01387     for (size_t i = 0; i < atoms->size(); i++) {
01388       eigenvec[i] -= eig_center;
01389     }
01390   }
01391 
01392   // eigenvec_invnorm2 is used when computing inverse gradients
01393   eigenvec_invnorm2 = 0.0;
01394   for (size_t ein = 0; ein < atoms->size(); ein++) {
01395     eigenvec_invnorm2 += eigenvec[ein].norm2();
01396   }
01397   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
01398 
01399   // Vector normalization overrides the default normalization for differenceVector
01400   bool normalize = false;
01401   get_keyval(conf, "normalizeVector", normalize, normalize);
01402 
01403   if (normalize) {
01404     cvm::log("Normalizing the vector so that |v| = 1.\n");
01405     for (size_t i = 0; i < atoms->size(); i++) {
01406       eigenvec[i] *= cvm::sqrt(eigenvec_invnorm2);
01407     }
01408     eigenvec_invnorm2 = 1.0;
01409   } else if (b_difference_vector) {
01410     cvm::log("Normalizing the vector so that the norm of the projection |v ⋅ (x_vec - x_ref)| = 1.\n");
01411     for (size_t i = 0; i < atoms->size(); i++) {
01412       eigenvec[i] *= eigenvec_invnorm2;
01413     }
01414     eigenvec_invnorm2 = 1.0/eigenvec_invnorm2;
01415   } else {
01416     cvm::log("The norm of the vector is |v| = "+
01417              cvm::to_str(1.0/cvm::sqrt(eigenvec_invnorm2))+".\n");
01418   }
01419 }
01420 
01421 
01422 void colvar::eigenvector::calc_value()
01423 {
01424   x.real_value = 0.0;
01425   for (size_t i = 0; i < atoms->size(); i++) {
01426     x.real_value += ((*atoms)[i].pos - ref_pos[i]) * eigenvec[i];
01427   }
01428 }
01429 
01430 
01431 void colvar::eigenvector::calc_gradients()
01432 {
01433   for (size_t ia = 0; ia < atoms->size(); ia++) {
01434     (*atoms)[ia].grad = eigenvec[ia];
01435   }
01436 }
01437 
01438 
01439 void colvar::eigenvector::apply_force(colvarvalue const &force)
01440 {
01441   if (!atoms->noforce)
01442     atoms->apply_colvar_force(force.real_value);
01443 }
01444 
01445 
01446 void colvar::eigenvector::calc_force_invgrads()
01447 {
01448   atoms->read_total_forces();
01449   ft.real_value = 0.0;
01450 
01451   for (size_t ia = 0; ia < atoms->size(); ia++) {
01452     ft.real_value += eigenvec_invnorm2 * (*atoms)[ia].grad *
01453       (*atoms)[ia].total_force;
01454   }
01455 }
01456 
01457 
01458 void colvar::eigenvector::calc_Jacobian_derivative()
01459 {
01460   // gradient of the rotation matrix
01461   cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
01462   cvm::quaternion &quat0 = atoms->rot.q;
01463 
01464   // gradients of products of 2 quaternion components
01465   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
01466 
01467   cvm::real sum = 0.0;
01468 
01469   for (size_t ia = 0; ia < atoms->size(); ia++) {
01470 
01471     // Gradient of optimal quaternion wrt current Cartesian position
01472     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
01473     // we can just transpose the derivatives of the direct matrix
01474     cvm::vector1d<cvm::rvector> &dq_1 = atoms->rot.dQ0_1[ia];
01475 
01476     g11 = 2.0 * quat0[1]*dq_1[1];
01477     g22 = 2.0 * quat0[2]*dq_1[2];
01478     g33 = 2.0 * quat0[3]*dq_1[3];
01479     g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
01480     g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
01481     g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
01482     g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
01483     g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
01484     g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
01485 
01486     // Gradient of the inverse rotation matrix wrt current Cartesian position
01487     // (transpose of the gradient of the direct rotation)
01488     grad_rot_mat[0][0] = -2.0 * (g22 + g33);
01489     grad_rot_mat[0][1] =  2.0 * (g12 + g03);
01490     grad_rot_mat[0][2] =  2.0 * (g13 - g02);
01491     grad_rot_mat[1][0] =  2.0 * (g12 - g03);
01492     grad_rot_mat[1][1] = -2.0 * (g11 + g33);
01493     grad_rot_mat[1][2] =  2.0 * (g01 + g23);
01494     grad_rot_mat[2][0] =  2.0 * (g02 + g13);
01495     grad_rot_mat[2][1] =  2.0 * (g23 - g01);
01496     grad_rot_mat[2][2] = -2.0 * (g11 + g22);
01497 
01498     for (size_t i = 0; i < 3; i++) {
01499       for (size_t j = 0; j < 3; j++) {
01500         sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
01501       }
01502     }
01503   }
01504 
01505   jd.real_value = sum * cvm::sqrt(eigenvec_invnorm2);
01506 }
01507 
01508 
01509 simple_scalar_dist_functions(eigenvector)
01510 
01511 
01512 
01513 colvar::cartesian::cartesian(std::string const &conf)
01514   : cvc(conf)
01515 {
01516   set_function_type("cartesian");
01517 
01518   atoms = parse_group(conf, "atoms");
01519 
01520   bool use_x, use_y, use_z;
01521   get_keyval(conf, "useX", use_x, true);
01522   get_keyval(conf, "useY", use_y, true);
01523   get_keyval(conf, "useZ", use_z, true);
01524 
01525   axes.clear();
01526   if (use_x) axes.push_back(0);
01527   if (use_y) axes.push_back(1);
01528   if (use_z) axes.push_back(2);
01529 
01530   if (axes.size() == 0) {
01531     cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
01532     return;
01533   }
01534 
01535   x.type(colvarvalue::type_vector);
01536   disable(f_cvc_explicit_gradient);
01537   // Don't try to access atoms if creation of the atom group failed
01538   if (atoms != NULL) x.vector1d_value.resize(atoms->size() * axes.size());
01539 }
01540 
01541 
01542 void colvar::cartesian::calc_value()
01543 {
01544   size_t const dim = axes.size();
01545   size_t ia, j;
01546   for (ia = 0; ia < atoms->size(); ia++) {
01547     for (j = 0; j < dim; j++) {
01548       x.vector1d_value[dim*ia + j] = (*atoms)[ia].pos[axes[j]];
01549     }
01550   }
01551 }
01552 
01553 
01554 void colvar::cartesian::calc_gradients()
01555 {
01556   // we're not using the "grad" member of each
01557   // atom object, because it only can represent the gradient of a
01558   // scalar colvar
01559 }
01560 
01561 
01562 void colvar::cartesian::apply_force(colvarvalue const &force)
01563 {
01564   size_t const dim = axes.size();
01565   size_t ia, j;
01566   if (!atoms->noforce) {
01567     cvm::rvector f;
01568     for (ia = 0; ia < atoms->size(); ia++) {
01569       for (j = 0; j < dim; j++) {
01570         f[axes[j]] = force.vector1d_value[dim*ia + j];
01571       }
01572       (*atoms)[ia].apply_force(f);
01573     }
01574   }
01575 }

Generated on Fri Oct 4 02:43:40 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002