colvarcomp_coordnums.C Source 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 "colvarmodule.h"
00011 #include "colvarparse.h"
00012 #include "colvaratoms.h"
00013 #include "colvarvalue.h"
00014 #include "colvar.h"
00015 #include "colvarcomp.h"
00016 
00017 
00018 
00019 template<int flags>
00020 cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
00021                                                cvm::rvector const &r0_vec,
00022                                                int en,
00023                                                int ed,
00024                                                cvm::atom &A1,
00025                                                cvm::atom &A2,
00026                                                bool **pairlist_elem,
00027                                                cvm::real pairlist_tol)
00028 {
00029   if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) {
00030     bool const within = **pairlist_elem;
00031     (*pairlist_elem)++;
00032     if (!within) {
00033       return 0.0;
00034     }
00035   }
00036 
00037   cvm::rvector const r0sq_vec(r0_vec.x*r0_vec.x,
00038                               r0_vec.y*r0_vec.y,
00039                               r0_vec.z*r0_vec.z);
00040 
00041   cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
00042 
00043   cvm::rvector const scal_diff(diff.x/((flags & ef_anisotropic) ?
00044                                        r0_vec.x : r0),
00045                                diff.y/((flags & ef_anisotropic) ?
00046                                        r0_vec.y : r0),
00047                                diff.z/((flags & ef_anisotropic) ?
00048                                        r0_vec.z : r0));
00049   cvm::real const l2 = scal_diff.norm2();
00050 
00051   // Assume en and ed are even integers, and avoid sqrt in the following
00052   int const en2 = en/2;
00053   int const ed2 = ed/2;
00054 
00055   cvm::real const xn = cvm::integer_power(l2, en2);
00056   cvm::real const xd = cvm::integer_power(l2, ed2);
00057   //The subtraction and division stretches the function back to the range of [0,1] from [pairlist_tol,1]
00058   cvm::real const func = (((1.0-xn)/(1.0-xd)) - pairlist_tol) / (1.0-pairlist_tol);
00059 
00060   if (flags & ef_rebuild_pairlist) {
00061     //Particles just outside of the cutoff also are considered if they come near.
00062     **pairlist_elem = (func > (-pairlist_tol * 0.5)) ? true : false;
00063     (*pairlist_elem)++;
00064   }
00065   //If the value is too small, we need to exclude it, rather than let it contribute to the sum or the gradients.
00066   if (func < 0)
00067     return 0;
00068 
00069   if (flags & ef_gradients) {
00070     //This is the old, completely correct expression for dFdl2:
00071     //cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
00072     //                                        func*ed2*(xd/l2))*(-1.0);
00073     //This can become:
00074     //cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2)*(1.0-xn)/(1.0-xn) -
00075     //                                        func*ed2*(xd/l2))*(-1.0);
00076     //Recognizing that func = (1.0-xn)/(1.0-xd), we can group together the "func" and get a version of dFdl2 that is 0
00077     //when func=0, which lets us skip this gradient calculation when func=0.
00078     cvm::real const dFdl2 = func * ((ed2*xd/((1.0-xd)*l2)) - (en2*xn/((1.0-xn)*l2)));
00079     cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
00080                                    r0*r0)) * diff.x,
00081                              (2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
00082                                    r0*r0)) * diff.y,
00083                              (2.0/((flags & ef_anisotropic) ? r0sq_vec.z :
00084                                    r0*r0)) * diff.z);
00085     A1.grad += (-1.0)*dFdl2*dl2dx;
00086     A2.grad +=        dFdl2*dl2dx;
00087   }
00088 
00089   return func;
00090 }
00091 
00092 
00093 colvar::coordnum::coordnum(std::string const &conf)
00094   : cvc(conf), b_anisotropic(false), pairlist(NULL)
00095 
00096 {
00097   set_function_type("coordNum");
00098   x.type(colvarvalue::type_scalar);
00099 
00100   colvarproxy *proxy = cvm::main()->proxy;
00101 
00102   group1 = parse_group(conf, "group1");
00103   group2 = parse_group(conf, "group2");
00104 
00105   if (group1 == NULL || group2 == NULL) {
00106     cvm::error("Error: failed to initialize atom groups.\n",
00107                 COLVARS_INPUT_ERROR);
00108     return;
00109   }
00110 
00111   if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
00112     cvm::error("Error: group1 and group2 share a common atom (number: " +
00113                cvm::to_str(atom_number) + ")\n", COLVARS_INPUT_ERROR);
00114     return;
00115   }
00116 
00117   if (group1->b_dummy) {
00118     cvm::error("Error: only group2 is allowed to be a dummy atom\n",
00119                COLVARS_INPUT_ERROR);
00120     return;
00121   }
00122 
00123   bool const b_isotropic = get_keyval(conf, "cutoff", r0,
00124                                       cvm::real(proxy->angstrom_to_internal(4.0)));
00125 
00126   if (get_keyval(conf, "cutoff3", r0_vec,
00127                  cvm::rvector(proxy->angstrom_to_internal(4.0),
00128                               proxy->angstrom_to_internal(4.0),
00129                               proxy->angstrom_to_internal(4.0)))) {
00130     if (b_isotropic) {
00131       cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
00132                  "at the same time.\n",
00133                  COLVARS_INPUT_ERROR);
00134       return;
00135     }
00136 
00137     b_anisotropic = true;
00138     // remove meaningless negative signs
00139     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
00140     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
00141     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
00142   }
00143 
00144   get_keyval(conf, "expNumer", en, 6);
00145   get_keyval(conf, "expDenom", ed, 12);
00146 
00147   if ( (en%2) || (ed%2) ) {
00148     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
00149                COLVARS_INPUT_ERROR);
00150   }
00151 
00152   if ( (en <= 0) || (ed <= 0) ) {
00153     cvm::error("Error: negative exponent(s) provided.\n",
00154                COLVARS_INPUT_ERROR);
00155   }
00156 
00157   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00158     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
00159   }
00160 
00161   get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
00162 
00163   get_keyval(conf, "tolerance", tolerance, 0.0);
00164   if (tolerance > 0) {
00165     cvm::main()->cite_feature("coordNum pairlist");
00166     get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
00167     if ( ! (pairlist_freq > 0) ) {
00168       cvm::error("Error: non-positive pairlistfrequency provided.\n",
00169                  COLVARS_INPUT_ERROR);
00170       return; // and do not allocate the pairlists below
00171     }
00172     if (b_group2_center_only) {
00173       pairlist = new bool[group1->size()];
00174     }
00175     else {
00176       pairlist = new bool[group1->size() * group2->size()];
00177     }
00178   }
00179 
00180   init_scalar_boundaries(0.0, b_group2_center_only ?
00181                          static_cast<cvm::real>(group1->size()) :
00182                          static_cast<cvm::real>(group1->size() *
00183                                                 group2->size()));
00184 }
00185 
00186 
00187 colvar::coordnum::~coordnum()
00188 {
00189   if (pairlist != NULL) {
00190     delete [] pairlist;
00191   }
00192 }
00193 
00194 
00195 template<int flags> void colvar::coordnum::main_loop(bool **pairlist_elem)
00196 {
00197   if (b_group2_center_only) {
00198     cvm::atom group2_com_atom;
00199     group2_com_atom.pos = group2->center_of_mass();
00200     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00201       x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
00202                                                 *ai1, group2_com_atom,
00203                                                 pairlist_elem,
00204                                                 tolerance);
00205     }
00206     if (b_group2_center_only) {
00207       group2->set_weighted_gradient(group2_com_atom.grad);
00208     }
00209   } else {
00210     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
00211       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
00212         x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
00213                                                   *ai1, *ai2,
00214                                                   pairlist_elem,
00215                                                   tolerance);
00216       }
00217     }
00218   }
00219 }
00220 
00221 
00222 template<int compute_flags> int colvar::coordnum::compute_coordnum()
00223 {
00224   bool const use_pairlist = (pairlist != NULL);
00225   bool const rebuild_pairlist = (pairlist != NULL) &&
00226     (cvm::step_relative() % pairlist_freq == 0);
00227 
00228   bool *pairlist_elem = use_pairlist ? pairlist : NULL;
00229 
00230   if (b_anisotropic) {
00231 
00232     if (use_pairlist) {
00233       if (rebuild_pairlist) {
00234         int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
00235           ef_rebuild_pairlist;
00236         main_loop<flags>(&pairlist_elem);
00237       } else {
00238         int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
00239         main_loop<flags>(&pairlist_elem);
00240       }
00241 
00242     } else {
00243 
00244       int const flags = compute_flags | ef_anisotropic;
00245       main_loop<flags>(NULL);
00246     }
00247 
00248   } else {
00249 
00250     if (use_pairlist) {
00251 
00252       if (rebuild_pairlist) {
00253         int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
00254         main_loop<flags>(&pairlist_elem);
00255       } else {
00256         int const flags = compute_flags | ef_use_pairlist;
00257         main_loop<flags>(&pairlist_elem);
00258       }
00259 
00260     } else {
00261 
00262       int const flags = compute_flags;
00263       main_loop<flags>(NULL);
00264     }
00265   }
00266 
00267   return COLVARS_OK;
00268 }
00269 
00270 
00271 void colvar::coordnum::calc_value()
00272 {
00273   x.real_value = 0.0;
00274   if (is_enabled(f_cvc_gradient)) {
00275     compute_coordnum<ef_gradients>();
00276   } else {
00277     compute_coordnum<ef_null>();
00278   }
00279 }
00280 
00281 
00282 void colvar::coordnum::calc_gradients()
00283 {
00284   // Gradients are computed by calc_value() if f_cvc_gradients is enabled
00285 }
00286 
00287 
00288 void colvar::coordnum::apply_force(colvarvalue const &force)
00289 {
00290   if (!group1->noforce)
00291     group1->apply_colvar_force(force.real_value);
00292 
00293   if (!group2->noforce)
00294     group2->apply_colvar_force(force.real_value);
00295 }
00296 
00297 
00298 simple_scalar_dist_functions(coordnum)
00299 
00300 
00301 
00302 // h_bond member functions
00303 
00304 colvar::h_bond::h_bond(std::string const &conf)
00305 : cvc(conf)
00306 {
00307   if (cvm::debug())
00308     cvm::log("Initializing h_bond object.\n");
00309 
00310   set_function_type("hBond");
00311   x.type(colvarvalue::type_scalar);
00312   init_scalar_boundaries(0.0, 1.0);
00313 
00314   colvarproxy *proxy = cvm::main()->proxy;
00315 
00316   int a_num = -1, d_num = -1;
00317   get_keyval(conf, "acceptor", a_num, a_num);
00318   get_keyval(conf, "donor",    d_num, a_num);
00319 
00320   if ( (a_num == -1) || (d_num == -1) ) {
00321     cvm::error("Error: either acceptor or donor undefined.\n");
00322     return;
00323   }
00324 
00325   cvm::atom acceptor = cvm::atom(a_num);
00326   cvm::atom donor    = cvm::atom(d_num);
00327   register_atom_group(new cvm::atom_group);
00328   atom_groups[0]->add_atom(acceptor);
00329   atom_groups[0]->add_atom(donor);
00330 
00331   get_keyval(conf, "cutoff",   r0, proxy->angstrom_to_internal(3.3));
00332   get_keyval(conf, "expNumer", en, 6);
00333   get_keyval(conf, "expDenom", ed, 8);
00334 
00335   if ( (en%2) || (ed%2) ) {
00336     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
00337                COLVARS_INPUT_ERROR);
00338   }
00339 
00340   if ( (en <= 0) || (ed <= 0) ) {
00341     cvm::error("Error: negative exponent(s) provided.\n",
00342                COLVARS_INPUT_ERROR);
00343   }
00344 
00345   if (cvm::debug())
00346     cvm::log("Done initializing h_bond object.\n");
00347 }
00348 
00349 
00350 colvar::h_bond::h_bond(cvm::atom const &acceptor,
00351                        cvm::atom const &donor,
00352                        cvm::real r0_i, int en_i, int ed_i)
00353   : r0(r0_i), en(en_i), ed(ed_i)
00354 {
00355   set_function_type("hBond");
00356   x.type(colvarvalue::type_scalar);
00357   init_scalar_boundaries(0.0, 1.0);
00358 
00359   register_atom_group(new cvm::atom_group);
00360   atom_groups[0]->add_atom(acceptor);
00361   atom_groups[0]->add_atom(donor);
00362 }
00363 
00364 
00365 void colvar::h_bond::calc_value()
00366 {
00367   int const flags = coordnum::ef_null;
00368   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
00369   x.real_value =
00370     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00371                                         (*atom_groups[0])[0],
00372                                         (*atom_groups[0])[1],
00373                                         NULL, 0.0);
00374 }
00375 
00376 
00377 void colvar::h_bond::calc_gradients()
00378 {
00379   int const flags = coordnum::ef_gradients;
00380   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
00381   coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00382                                       (*atom_groups[0])[0],
00383                                       (*atom_groups[0])[1],
00384                                       NULL, 0.0);
00385 }
00386 
00387 
00388 void colvar::h_bond::apply_force(colvarvalue const &force)
00389 {
00390   (atom_groups[0])->apply_colvar_force(force);
00391 }
00392 
00393 
00394 simple_scalar_dist_functions(h_bond)
00395 
00396 
00397 
00398 colvar::selfcoordnum::selfcoordnum(std::string const &conf)
00399   : cvc(conf), pairlist(NULL)
00400 {
00401   set_function_type("selfCoordNum");
00402   x.type(colvarvalue::type_scalar);
00403 
00404   colvarproxy *proxy = cvm::main()->proxy;
00405 
00406   group1 = parse_group(conf, "group1");
00407 
00408   get_keyval(conf, "cutoff", r0, cvm::real(proxy->angstrom_to_internal(4.0)));
00409   get_keyval(conf, "expNumer", en, 6);
00410   get_keyval(conf, "expDenom", ed, 12);
00411 
00412 
00413   if ( (en%2) || (ed%2) ) {
00414     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
00415                COLVARS_INPUT_ERROR);
00416   }
00417 
00418   if ( (en <= 0) || (ed <= 0) ) {
00419     cvm::error("Error: negative exponent(s) provided.\n",
00420                COLVARS_INPUT_ERROR);
00421   }
00422 
00423   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00424     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
00425   }
00426 
00427   get_keyval(conf, "tolerance", tolerance, 0.0);
00428   if (tolerance > 0) {
00429     get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
00430     if ( ! (pairlist_freq > 0) ) {
00431       cvm::error("Error: non-positive pairlistfrequency provided.\n",
00432                  COLVARS_INPUT_ERROR);
00433       return;
00434     }
00435     pairlist = new bool[(group1->size()-1) * (group1->size()-1)];
00436   }
00437 
00438   init_scalar_boundaries(0.0, static_cast<cvm::real>((group1->size()-1) *
00439                                                      (group1->size()-1)));
00440 }
00441 
00442 
00443 colvar::selfcoordnum::~selfcoordnum()
00444 {
00445   if (pairlist != NULL) {
00446     delete [] pairlist;
00447   }
00448 }
00449 
00450 
00451 template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
00452 {
00453   cvm::rvector const r0_vec(0.0); // TODO enable the flag?
00454 
00455   bool const use_pairlist = (pairlist != NULL);
00456   bool const rebuild_pairlist = (pairlist != NULL) &&
00457     (cvm::step_relative() % pairlist_freq == 0);
00458 
00459   bool *pairlist_elem = use_pairlist ? pairlist : NULL;
00460   size_t i = 0, j = 0;
00461   size_t const n = group1->size();
00462 
00463   // Always isotropic (TODO: enable the ellipsoid?)
00464 
00465   if (use_pairlist) {
00466 
00467     if (rebuild_pairlist) {
00468       int const flags = compute_flags | coordnum::ef_use_pairlist |
00469         coordnum::ef_rebuild_pairlist;
00470       for (i = 0; i < n - 1; i++) {
00471         for (j = i + 1; j < n; j++) {
00472           x.real_value +=
00473             coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00474                                                 (*group1)[i],
00475                                                 (*group1)[j],
00476                                                 &pairlist_elem,
00477                                                 tolerance);
00478         }
00479       }
00480     } else {
00481       int const flags = compute_flags | coordnum::ef_use_pairlist;
00482       for (i = 0; i < n - 1; i++) {
00483         for (j = i + 1; j < n; j++) {
00484           x.real_value +=
00485             coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00486                                                 (*group1)[i],
00487                                                 (*group1)[j],
00488                                                 &pairlist_elem,
00489                                                 tolerance);
00490         }
00491       }
00492     }
00493 
00494   } else { // if (use_pairlist) {
00495 
00496     int const flags = compute_flags | coordnum::ef_null;
00497     for (i = 0; i < n - 1; i++) {
00498       for (j = i + 1; j < n; j++) {
00499         x.real_value +=
00500           coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00501                                               (*group1)[i],
00502                                               (*group1)[j],
00503                                               &pairlist_elem,
00504                                               tolerance);
00505       }
00506     }
00507   }
00508 
00509   return COLVARS_OK;
00510 }
00511 
00512 
00513 void colvar::selfcoordnum::calc_value()
00514 {
00515   x.real_value = 0.0;
00516   if (is_enabled(f_cvc_gradient)) {
00517     compute_selfcoordnum<coordnum::ef_gradients>();
00518   } else {
00519     compute_selfcoordnum<coordnum::ef_null>();
00520   }
00521 }
00522 
00523 
00524 void colvar::selfcoordnum::calc_gradients()
00525 {
00526   // Gradients are computed by calc_value() if f_cvc_gradients is enabled
00527 }
00528 
00529 
00530 void colvar::selfcoordnum::apply_force(colvarvalue const &force)
00531 {
00532   if (!group1->noforce) {
00533     group1->apply_colvar_force(force.real_value);
00534   }
00535 }
00536 
00537 
00538 simple_scalar_dist_functions(selfcoordnum)
00539 
00540 
00541 
00542 // groupcoordnum member functions
00543 colvar::groupcoordnum::groupcoordnum(std::string const &conf)
00544   : distance(conf), b_anisotropic(false)
00545 {
00546   set_function_type("groupCoord");
00547   x.type(colvarvalue::type_scalar);
00548   init_scalar_boundaries(0.0, 1.0);
00549 
00550   colvarproxy *proxy = cvm::main()->proxy;
00551 
00552   // group1 and group2 are already initialized by distance()
00553   if (group1->b_dummy || group2->b_dummy) {
00554     cvm::error("Error: neither group can be a dummy atom\n");
00555     return;
00556   }
00557 
00558   bool const b_scale = get_keyval(conf, "cutoff", r0,
00559                                   cvm::real(proxy->angstrom_to_internal(4.0)));
00560 
00561   if (get_keyval(conf, "cutoff3", r0_vec,
00562                  cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
00563 
00564     if (b_scale) {
00565       cvm::error("Error: cannot specify \"scale\" and "
00566                  "\"scale3\" at the same time.\n");
00567       return;
00568     }
00569     b_anisotropic = true;
00570     // remove meaningless negative signs
00571     if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
00572     if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
00573     if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
00574   }
00575 
00576   get_keyval(conf, "expNumer", en, 6);
00577   get_keyval(conf, "expDenom", ed, 12);
00578 
00579   if ( (en%2) || (ed%2) ) {
00580     cvm::error("Error: odd exponent(s) provided, can only use even ones.\n",
00581                COLVARS_INPUT_ERROR);
00582   }
00583 
00584   if ( (en <= 0) || (ed <= 0) ) {
00585     cvm::error("Error: negative exponent(s) provided.\n",
00586                COLVARS_INPUT_ERROR);
00587   }
00588 
00589   if (!is_enabled(f_cvc_pbc_minimum_image)) {
00590     cvm::log("Warning: only minimum-image distances are used by this variable.\n");
00591   }
00592 
00593 }
00594 
00595 
00596 void colvar::groupcoordnum::calc_value()
00597 {
00598   // create fake atoms to hold the com coordinates
00599   cvm::atom group1_com_atom;
00600   cvm::atom group2_com_atom;
00601   group1_com_atom.pos = group1->center_of_mass();
00602   group2_com_atom.pos = group2->center_of_mass();
00603   if (b_anisotropic) {
00604     int const flags = coordnum::ef_anisotropic;
00605     x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00606                                                        group1_com_atom,
00607                                                        group2_com_atom,
00608                                                        NULL, 0.0);
00609   } else {
00610     int const flags = coordnum::ef_null;
00611     x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00612                                                        group1_com_atom,
00613                                                        group2_com_atom,
00614                                                        NULL, 0.0);
00615   }
00616 }
00617 
00618 
00619 void colvar::groupcoordnum::calc_gradients()
00620 {
00621   cvm::atom group1_com_atom;
00622   cvm::atom group2_com_atom;
00623   group1_com_atom.pos = group1->center_of_mass();
00624   group2_com_atom.pos = group2->center_of_mass();
00625 
00626   if (b_anisotropic) {
00627     int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
00628     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00629                                         group1_com_atom,
00630                                         group2_com_atom,
00631                                         NULL, 0.0);
00632   } else {
00633     int const flags = coordnum::ef_gradients;
00634     coordnum::switching_function<flags>(r0, r0_vec, en, ed,
00635                                         group1_com_atom,
00636                                         group2_com_atom,
00637                                         NULL, 0.0);
00638   }
00639 
00640   group1->set_weighted_gradient(group1_com_atom.grad);
00641   group2->set_weighted_gradient(group2_com_atom.grad);
00642 }
00643 
00644 
00645 void colvar::groupcoordnum::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(groupcoordnum)