Difference for src/colvarcomp_distances.C from version 1.15 to 1.16

version 1.15version 1.16
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include <cmath> #include <cmath>
  
Line 10
Line 10
  
  
  
 // "twogroup" flag defaults to true; set to false by selfCoordNum colvar::distance::distance(std::string const &conf)
 // (only distance-derived component based on only one group) 
  
 colvar::distance::distance(std::string const &conf, bool twogroups) 
   : cvc(conf)   : cvc(conf)
 { {
   function_type = "distance";   function_type = "distance";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
  
    group1 = parse_group(conf, "group1");
    group2 = parse_group(conf, "group2");
  
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }   }
   if (twogroups && get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { 
     cvm::log("Computing system force on group 1 only");   init_total_force_params(conf);
   } 
   parse_group(conf, "group1", group1); 
   atom_groups.push_back(&group1); 
   if (twogroups) { 
     parse_group(conf, "group2", group2); 
     atom_groups.push_back(&group2); 
   } 
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 39
Line 35
   : cvc()   : cvc()
 { {
   function_type = "distance";   function_type = "distance";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   b_1site_force = false;   provide(f_cvc_com_based);
   b_no_PBC = false;   b_no_PBC = false;
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
  
 void colvar::distance::calc_value() void colvar::distance::calc_value()
 { {
   if (b_no_PBC) {   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();     dist_v = group2->center_of_mass() - group1->center_of_mass();
   } else {   } else {
     dist_v = cvm::position_distance(group1.center_of_mass(),     dist_v = cvm::position_distance(group1->center_of_mass(),
                                      group2.center_of_mass());                                      group2->center_of_mass());
   }   }
   x.real_value = dist_v.norm();   x.real_value = dist_v.norm();
 } }
  
  
 void colvar::distance::calc_gradients() void colvar::distance::calc_gradients()
 { {
   cvm::rvector const u = dist_v.unit();   cvm::rvector const u = dist_v.unit();
   group1.set_weighted_gradient(-1.0 * u);   group1->set_weighted_gradient(-1.0 * u);
   group2.set_weighted_gradient(       u);   group2->set_weighted_gradient(       u);
 } }
  
  
 void colvar::distance::calc_force_invgrads() void colvar::distance::calc_force_invgrads()
 { {
   group1.read_system_forces();   group1->read_total_forces();
   if ( b_1site_force ) {   if (is_enabled(f_cvc_one_site_total_force)) {
     ft.real_value = -1.0 * (group1.system_force() * dist_v.unit());     ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
   } else {   } else {
     group2.read_system_forces();     group2->read_total_forces();
     ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit());     ft.real_value = 0.5 * ((group2->total_force() - group1->total_force()) * dist_v.unit());
   }   }
 } }
  
  
 void colvar::distance::calc_Jacobian_derivative() void colvar::distance::calc_Jacobian_derivative()
 { {
   jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;   jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
 } }
  
  
 void colvar::distance::apply_force(colvarvalue const &force) void colvar::distance::apply_force(colvarvalue const &force)
 { {
   if (!group1.noforce)   if (!group1->noforce)
     group1.apply_colvar_force(force.real_value);     group1->apply_colvar_force(force.real_value);
  
   if (!group2.noforce)   if (!group2->noforce)
     group2.apply_colvar_force(force.real_value);     group2->apply_colvar_force(force.real_value);
 } }
  
  
Line 95
Line 96
   : distance(conf)   : distance(conf)
 { {
   function_type = "distance_vec";   function_type = "distance_vec";
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);   x.type(colvarvalue::type_3vector);
 } }
  
  
 colvar::distance_vec::distance_vec() colvar::distance_vec::distance_vec()
   : distance()   : distance()
 { {
   function_type = "distance_vec";   function_type = "distance_vec";
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);   x.type(colvarvalue::type_3vector);
 } }
  
  
 void colvar::distance_vec::calc_value() void colvar::distance_vec::calc_value()
 { {
   if (b_no_PBC) {   if (b_no_PBC) {
     x.rvector_value = group2.center_of_mass() - group1.center_of_mass();     x.rvector_value = group2->center_of_mass() - group1->center_of_mass();
   } else {   } else {
     x.rvector_value = cvm::position_distance(group1.center_of_mass(),     x.rvector_value = cvm::position_distance(group1->center_of_mass(),
                                               group2.center_of_mass());                                               group2->center_of_mass());
   }   }
 } }
  
  
 void colvar::distance_vec::calc_gradients() void colvar::distance_vec::calc_gradients()
 { {
   // gradients are not stored: a 3x3 matrix for each atom would be   // gradients are not stored: a 3x3 matrix for each atom would be
   // needed to store just the identity matrix   // needed to store just the identity matrix
 } }
  
  
 void colvar::distance_vec::apply_force(colvarvalue const &force) void colvar::distance_vec::apply_force(colvarvalue const &force)
 { {
   if (!group1.noforce)   if (!group1->noforce)
     group1.apply_force(-1.0 * force.rvector_value);     group1->apply_force(-1.0 * force.rvector_value);
  
   if (!group2.noforce)   if (!group2->noforce)
     group2.apply_force(       force.rvector_value);     group2->apply_force(       force.rvector_value);
 } }
  
  
Line 136
Line 143
   : cvc(conf)   : cvc(conf)
 { {
   function_type = "distance_z";   function_type = "distance_z";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
   // TODO detect PBC from MD engine (in simple cases)   // TODO detect PBC from MD engine (in simple cases)
Line 151
Line 159
     return;     return;
   }   }
  
   parse_group(conf, "main", main);   main = parse_group(conf, "main");
   parse_group(conf, "ref", ref1);   ref1 = parse_group(conf, "ref");
   atom_groups.push_back(&main); 
   atom_groups.push_back(&ref1); 
   // this group is optional   // this group is optional
   parse_group(conf, "ref2", ref2, true);   ref2 = parse_group(conf, "ref2", true);
  
   if (ref2.size()) {   if (ref2 && ref2->size()) {
     atom_groups.push_back(&ref2); 
     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
     fixed_axis = false;     fixed_axis = false;
     if (key_lookup(conf, "axis"))     if (key_lookup(conf, "axis"))
Line 181
Line 186
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }   }
   if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { 
     cvm::log("Computing system force on group \"main\" only");   init_total_force_params(conf);
   } 
 } }
  
 colvar::distance_z::distance_z() colvar::distance_z::distance_z()
 { {
   function_type = "distance_z";   function_type = "distance_z";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 198
Line 204
 { {
   if (fixed_axis) {   if (fixed_axis) {
     if (b_no_PBC) {     if (b_no_PBC) {
       dist_v = main.center_of_mass() - ref1.center_of_mass();       dist_v = main->center_of_mass() - ref1->center_of_mass();
     } else {     } else {
       dist_v = cvm::position_distance(ref1.center_of_mass(),       dist_v = cvm::position_distance(ref1->center_of_mass(),
                                        main.center_of_mass());                                        main->center_of_mass());
     }     }
   } else {   } else {
  
     if (b_no_PBC) {     if (b_no_PBC) {
       dist_v = main.center_of_mass() -       dist_v = main->center_of_mass() -
                (0.5 * (ref1.center_of_mass() + ref2.center_of_mass()));                (0.5 * (ref1->center_of_mass() + ref2->center_of_mass()));
       axis = ref2.center_of_mass() - ref1.center_of_mass();       axis = ref2->center_of_mass() - ref1->center_of_mass();
     } else {     } else {
       dist_v = cvm::position_distance(0.5 * (ref1.center_of_mass() +       dist_v = cvm::position_distance(0.5 * (ref1->center_of_mass() +
                ref2.center_of_mass()), main.center_of_mass());                ref2->center_of_mass()), main->center_of_mass());
       axis = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass());       axis = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass());
     }     }
     axis_norm = axis.norm();     axis_norm = axis.norm();
     axis = axis.unit();     axis = axis.unit();
Line 223
Line 229
  
 void colvar::distance_z::calc_gradients() void colvar::distance_z::calc_gradients()
 { {
   main.set_weighted_gradient( axis );   main->set_weighted_gradient( axis );
  
   if (fixed_axis) {   if (fixed_axis) {
     ref1.set_weighted_gradient(-1.0 * axis);     ref1->set_weighted_gradient(-1.0 * axis);
   } else {   } else {
     if (b_no_PBC) {     if (b_no_PBC) {
       ref1.set_weighted_gradient( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() -       ref1->set_weighted_gradient( 1.0 / axis_norm * (main->center_of_mass() - ref2->center_of_mass() -
                                    x.real_value * axis ));                                    x.real_value * axis ));
       ref2.set_weighted_gradient( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() +       ref2->set_weighted_gradient( 1.0 / axis_norm * (ref1->center_of_mass() - main->center_of_mass() +
                                    x.real_value * axis ));                                    x.real_value * axis ));
     } else {     } else {
       ref1.set_weighted_gradient( 1.0 / axis_norm * (       ref1->set_weighted_gradient( 1.0 / axis_norm * (
         cvm::position_distance(ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis ));         cvm::position_distance(ref2->center_of_mass(), main->center_of_mass()) - x.real_value * axis ));
       ref2.set_weighted_gradient( 1.0 / axis_norm * (       ref2->set_weighted_gradient( 1.0 / axis_norm * (
         cvm::position_distance(main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis ));         cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis ));
     } 
   }   }
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients for group main:\n"); 
     debug_gradients(main); 
     cvm::log("Debugging gradients for group ref1:\n"); 
     debug_gradients(ref1); 
     cvm::log("Debugging gradients for group ref2:\n"); 
     debug_gradients(ref2); 
   }   }
 } }
  
 void colvar::distance_z::calc_force_invgrads() void colvar::distance_z::calc_force_invgrads()
 { {
   main.read_system_forces();   main->read_total_forces();
  
   if (fixed_axis && !b_1site_force) {   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1.read_system_forces();     ref1->read_total_forces();
     ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis);     ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
   } else {   } else {
     ft.real_value = main.system_force() * axis;     ft.real_value = main->total_force() * axis;
   }   }
 } }
  
Line 270
Line 267
  
 void colvar::distance_z::apply_force(colvarvalue const &force) void colvar::distance_z::apply_force(colvarvalue const &force)
 { {
   if (!ref1.noforce)   if (!ref1->noforce)
     ref1.apply_colvar_force(force.real_value);     ref1->apply_colvar_force(force.real_value);
  
   if (ref2.size() && !ref2.noforce)   if (ref2 && ref2->size() && !ref2->noforce)
     ref2.apply_colvar_force(force.real_value);     ref2->apply_colvar_force(force.real_value);
  
   if (!main.noforce)   if (!main->noforce)
     main.apply_colvar_force(force.real_value);     main->apply_colvar_force(force.real_value);
 } }
  
  
Line 286
Line 283
   : distance_z(conf)   : distance_z(conf)
 { {
   function_type = "distance_xy";   function_type = "distance_xy";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 295
Line 293
   : distance_z()   : distance_z()
 { {
   function_type = "distance_xy";   function_type = "distance_xy";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
 void colvar::distance_xy::calc_value() void colvar::distance_xy::calc_value()
 { {
   if (b_no_PBC) {   if (b_no_PBC) {
     dist_v = main.center_of_mass() - ref1.center_of_mass();     dist_v = main->center_of_mass() - ref1->center_of_mass();
   } else {   } else {
     dist_v = cvm::position_distance(ref1.center_of_mass(),     dist_v = cvm::position_distance(ref1->center_of_mass(),
                                      main.center_of_mass());                                      main->center_of_mass());
   }   }
   if (!fixed_axis) {   if (!fixed_axis) {
     if (b_no_PBC) {     if (b_no_PBC) {
       v12 = ref2.center_of_mass() - ref1.center_of_mass();       v12 = ref2->center_of_mass() - ref1->center_of_mass();
     } else {     } else {
       v12 = cvm::position_distance(ref1.center_of_mass(), ref2.center_of_mass());       v12 = cvm::position_distance(ref1->center_of_mass(), ref2->center_of_mass());
     }     }
     axis_norm = v12.norm();     axis_norm = v12.norm();
     axis = v12.unit();     axis = v12.unit();
Line 333
Line 332
   x_inv = 1.0 / x.real_value;   x_inv = 1.0 / x.real_value;
  
   if (fixed_axis) {   if (fixed_axis) {
     ref1.set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);     ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho);
     main.set_weighted_gradient(       x_inv * dist_v_ortho);     main->set_weighted_gradient(       x_inv * dist_v_ortho);
   } else {   } else {
     if (b_no_PBC) {     if (b_no_PBC) {
       v13 = main.center_of_mass() - ref1.center_of_mass();       v13 = main->center_of_mass() - ref1->center_of_mass();
     } else {     } else {
       v13 = cvm::position_distance(ref1.center_of_mass(), main.center_of_mass());       v13 = cvm::position_distance(ref1->center_of_mass(), main->center_of_mass());
     }     }
     A = (dist_v * axis) / axis_norm;     A = (dist_v * axis) / axis_norm;
  
     ref1.set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho);     ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho);
     ref2.set_weighted_gradient( -A        * x_inv * dist_v_ortho);     ref2->set_weighted_gradient( -A        * x_inv * dist_v_ortho);
     main.set_weighted_gradient(      1.0  * x_inv * dist_v_ortho);     main->set_weighted_gradient(      1.0  * x_inv * dist_v_ortho);
   }   }
 } }
  
 void colvar::distance_xy::calc_force_invgrads() void colvar::distance_xy::calc_force_invgrads()
 { {
   main.read_system_forces();   main->read_total_forces();
  
   if (fixed_axis && !b_1site_force) {   if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1.read_system_forces();     ref1->read_total_forces();
     ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho);     ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
   } else {   } else {
     ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho;     ft.real_value = 1.0 / x.real_value * main->total_force() * dist_v_ortho;
   }   }
 } }
  
Line 368
Line 367
  
 void colvar::distance_xy::apply_force(colvarvalue const &force) void colvar::distance_xy::apply_force(colvarvalue const &force)
 { {
   if (!ref1.noforce)   if (!ref1->noforce)
     ref1.apply_colvar_force(force.real_value);     ref1->apply_colvar_force(force.real_value);
  
   if (ref2.size() && !ref2.noforce)   if (ref2 && ref2->size() && !ref2->noforce)
     ref2.apply_colvar_force(force.real_value);     ref2->apply_colvar_force(force.real_value);
  
   if (!main.noforce)   if (!main->noforce)
     main.apply_colvar_force(force.real_value);     main->apply_colvar_force(force.real_value);
 } }
  
  
Line 384
Line 383
   : distance(conf)   : distance(conf)
 { {
   function_type = "distance_dir";   function_type = "distance_dir";
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);   x.type(colvarvalue::type_unit3vector);
 } }
  
Line 392
Line 392
   : distance()   : distance()
 { {
   function_type = "distance_dir";   function_type = "distance_dir";
    provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);   x.type(colvarvalue::type_unit3vector);
 } }
  
Line 399
Line 400
 void colvar::distance_dir::calc_value() void colvar::distance_dir::calc_value()
 { {
   if (b_no_PBC) {   if (b_no_PBC) {
     dist_v = group2.center_of_mass() - group1.center_of_mass();     dist_v = group2->center_of_mass() - group1->center_of_mass();
   } else {   } else {
     dist_v = cvm::position_distance(group1.center_of_mass(),     dist_v = cvm::position_distance(group1->center_of_mass(),
                                      group2.center_of_mass());                                      group2->center_of_mass());
   }   }
   x.rvector_value = dist_v.unit();   x.rvector_value = dist_v.unit();
 } }
Line 413
Line 414
   // gradients are computed on the fly within apply_force()   // gradients are computed on the fly within apply_force()
   // Note: could be a problem if a future bias relies on gradient   // Note: could be a problem if a future bias relies on gradient
   // calculations...   // calculations...
    // TODO in new deps system: remove dependency of biasing force to gradient?
    // That way we could tell apart an explicit gradient dependency
 } }
  
  
Line 422
Line 425
   cvm::real const iprod = force.rvector_value * x.rvector_value;   cvm::real const iprod = force.rvector_value * x.rvector_value;
   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;   cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
  
   if (!group1.noforce)   if (!group1->noforce)
     group1.apply_force(-1.0 * force_tang);     group1->apply_force(-1.0 * force_tang);
  
   if (!group2.noforce)   if (!group2->noforce)
     group2.apply_force(       force_tang);     group2->apply_force(       force_tang);
 } }
  
  
Line 445
Line 448
     return;     return;
   }   }
  
   for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
     for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {     for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
       if (ai1->id == ai2->id) {       if (ai1->id == ai2->id) {
         cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n");         cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n");
         return;         return;
Line 454
Line 457
     }     }
   }   }
  
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 463
Line 464
 { {
   function_type = "distance_inv";   function_type = "distance_inv";
   exponent = 6;   exponent = 6;
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
   b_1site_force = false; 
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 473
Line 471
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   if (b_no_PBC) {   if (b_no_PBC) {
     for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
         cvm::rvector const dv = ai2->pos - ai1->pos;         cvm::rvector const dv = ai2->pos - ai1->pos;
         cvm::real const d2 = dv.norm2();         cvm::real const d2 = dv.norm2();
         cvm::real dinv = 1.0;         cvm::real dinv = 1.0;
Line 487
Line 485
       }       }
     }     }
   } else {   } else {
     for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {     for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {       for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
         cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos);         cvm::rvector const dv = cvm::position_distance(ai1->pos, ai2->pos);
         cvm::real const d2 = dv.norm2();         cvm::real const d2 = dv.norm2();
         cvm::real dinv = 1.0;         cvm::real dinv = 1.0;
Line 502
Line 500
     }     }
   }   }
  
   x.real_value *= 1.0 / cvm::real(group1.size() * group2.size());   x.real_value *= 1.0 / cvm::real(group1->size() * group2->size());
   x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));   x.real_value = std::pow(x.real_value, -1.0/(cvm::real(exponent)));
 } }
  
 void colvar::distance_inv::calc_gradients() void colvar::distance_inv::calc_gradients()
 { {
   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1.size() * group2.size());   cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * std::pow(x.real_value, exponent+1) / cvm::real(group1->size() * group2->size());
   for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {   for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
     ai1->grad *= dxdsum;     ai1->grad *= dxdsum;
   }   }
   for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {   for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
     ai2->grad *= dxdsum;     ai2->grad *= dxdsum;
   }   }
 } }
  
 void colvar::distance_inv::apply_force(colvarvalue const &force) void colvar::distance_inv::apply_force(colvarvalue const &force)
 { {
   if (!group1.noforce)   if (!group1->noforce)
     group1.apply_colvar_force(force.real_value);     group1->apply_colvar_force(force.real_value);
  
   if (!group2.noforce)   if (!group2->noforce)
     group2.apply_colvar_force(force.real_value);     group2->apply_colvar_force(force.real_value);
 } }
  
  
Line 532
Line 530
   : cvc(conf)   : cvc(conf)
 { {
   function_type = "distance_pairs";   function_type = "distance_pairs";
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
  
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }   }
  
   parse_group(conf, "group1", group1);   group1 = parse_group(conf, "group1");
   atom_groups.push_back(&group1);   group2 = parse_group(conf, "group2");
  
   parse_group(conf, "group2", group2); 
   atom_groups.push_back(&group2); 
  
   x.type(colvarvalue::type_vector);   x.type(colvarvalue::type_vector);
   x.vector1d_value.resize(group1.size() * group2.size());   x.vector1d_value.resize(group1->size() * group2->size());
 } }
  
  
 colvar::distance_pairs::distance_pairs() colvar::distance_pairs::distance_pairs()
 { {
   function_type = "distance_pairs";   function_type = "distance_pairs";
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
   x.type(colvarvalue::type_vector);   x.type(colvarvalue::type_vector);
 } }
  
  
 void colvar::distance_pairs::calc_value() void colvar::distance_pairs::calc_value()
 { {
   x.vector1d_value.resize(group1.size() * group2.size());   x.vector1d_value.resize(group1->size() * group2->size());
  
   if (b_no_PBC) {   if (b_no_PBC) {
     size_t i1, i2;     size_t i1, i2;
     for (i1 = 0; i1 < group1.size(); i1++) {     for (i1 = 0; i1 < group1->size(); i1++) {
       for (i2 = 0; i2 < group2.size(); i2++) {       for (i2 = 0; i2 < group2->size(); i2++) {
         cvm::rvector const dv = group2[i2].pos - group1[i1].pos;         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
         cvm::real const d = dv.norm();         cvm::real const d = dv.norm();
         x.vector1d_value[i1*group2.size() + i2] = d;         x.vector1d_value[i1*group2->size() + i2] = d;
         group1[i1].grad = -1.0 * dv.unit();         (*group1)[i1].grad = -1.0 * dv.unit();
         group2[i2].grad =  dv.unit();         (*group2)[i2].grad =  dv.unit();
       }       }
     }     }
   } else {   } else {
     size_t i1, i2;     size_t i1, i2;
     for (i1 = 0; i1 < group1.size(); i1++) {     for (i1 = 0; i1 < group1->size(); i1++) {
       for (i2 = 0; i2 < group2.size(); i2++) {       for (i2 = 0; i2 < group2->size(); i2++) {
         cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos);         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos);
         cvm::real const d = dv.norm();         cvm::real const d = dv.norm();
         x.vector1d_value[i1*group2.size() + i2] = d;         x.vector1d_value[i1*group2->size() + i2] = d;
         group1[i1].grad = -1.0 * dv.unit();         (*group1)[i1].grad = -1.0 * dv.unit();
         group2[i2].grad =  dv.unit();         (*group2)[i2].grad =  dv.unit();
       }       }
     }     }
   }   }
Line 591
Line 582
 void colvar::distance_pairs::calc_gradients() void colvar::distance_pairs::calc_gradients()
 { {
   // will be calculated on the fly in apply_force()   // will be calculated on the fly in apply_force()
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(group1); 
   } 
 } }
  
 void colvar::distance_pairs::apply_force(colvarvalue const &force) void colvar::distance_pairs::apply_force(colvarvalue const &force)
 { {
   if (b_no_PBC) {   if (b_no_PBC) {
     size_t i1, i2;     size_t i1, i2;
     for (i1 = 0; i1 < group1.size(); i1++) {     for (i1 = 0; i1 < group1->size(); i1++) {
       for (i2 = 0; i2 < group2.size(); i2++) {       for (i2 = 0; i2 < group2->size(); i2++) {
         cvm::rvector const dv = group2[i2].pos - group1[i1].pos;         cvm::rvector const dv = (*group2)[i2].pos - (*group1)[i1].pos;
         group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit());         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
         group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit());         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
       }       }
     }     }
   } else {   } else {
     size_t i1, i2;     size_t i1, i2;
     for (i1 = 0; i1 < group1.size(); i1++) {     for (i1 = 0; i1 < group1->size(); i1++) {
       for (i2 = 0; i2 < group2.size(); i2++) {       for (i2 = 0; i2 < group2->size(); i2++) {
         cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos);         cvm::rvector const dv = cvm::position_distance((*group1)[i1].pos, (*group2)[i2].pos);
         group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit());         (*group1)[i1].apply_force(force[i1*group2->size() + i2] * (-1.0) * dv.unit());
         group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit());         (*group2)[i2].apply_force(force[i1*group2->size() + i2] * dv.unit());
       }       }
     }     }
   }   }
Line 625
Line 612
   : cvc(conf)   : cvc(conf)
 { {
   function_type = "gyration";   function_type = "gyration";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   parse_group(conf, "atoms", atoms);   atoms = parse_group(conf, "atoms");
   atom_groups.push_back(&atoms); 
  
   if (atoms.b_user_defined_fit) {   if (atoms->b_user_defined_fit) {
     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
   } else {   } else {
     atoms.b_center = true;     atoms->b_center = true;
     atoms.ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));     atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0));
   }   }
  
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
Line 644
Line 630
 colvar::gyration::gyration() colvar::gyration::gyration()
 { {
   function_type = "gyration";   function_type = "gyration";
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 653
Line 639
 void colvar::gyration::calc_value() void colvar::gyration::calc_value()
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     x.real_value += (ai->pos).norm2();     x.real_value += (ai->pos).norm2();
   }   }
   x.real_value = std::sqrt(x.real_value / cvm::real(atoms.size()));   x.real_value = std::sqrt(x.real_value / cvm::real(atoms->size()));
 } }
  
  
 void colvar::gyration::calc_gradients() void colvar::gyration::calc_gradients()
 { {
   cvm::real const drdx = 1.0/(cvm::real(atoms.size()) * x.real_value);   cvm::real const drdx = 1.0/(cvm::real(atoms->size()) * x.real_value);
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = drdx * ai->pos;     ai->grad = drdx * ai->pos;
   }   }
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(atoms); 
   } 
 } }
  
  
 void colvar::gyration::calc_force_invgrads() void colvar::gyration::calc_force_invgrads()
 { {
   atoms.read_system_forces();   atoms->read_total_forces();
  
   cvm::real const dxdr = 1.0/x.real_value;   cvm::real const dxdr = 1.0/x.real_value;
   ft.real_value = 0.0;   ft.real_value = 0.0;
  
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ft.real_value += dxdr * ai->pos * ai->system_force;     ft.real_value += dxdr * ai->pos * ai->total_force;
   }   }
 } }
  
  
 void colvar::gyration::calc_Jacobian_derivative() void colvar::gyration::calc_Jacobian_derivative()
 { {
   jd = x.real_value ? (3.0 * cvm::real(atoms.size()) - 4.0) / x.real_value : 0.0;   jd = x.real_value ? (3.0 * cvm::real(atoms->size()) - 4.0) / x.real_value : 0.0;
 } }
  
  
 void colvar::gyration::apply_force(colvarvalue const &force) void colvar::gyration::apply_force(colvarvalue const &force)
 { {
   if (!atoms.noforce)   if (!atoms->noforce)
     atoms.apply_colvar_force(force.real_value);     atoms->apply_colvar_force(force.real_value);
 } }
  
  
Line 705
Line 686
   : gyration(conf)   : gyration(conf)
 { {
   function_type = "inertia";   function_type = "inertia";
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
 } }
  
Line 721
Line 700
 void colvar::inertia::calc_value() void colvar::inertia::calc_value()
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     x.real_value += (ai->pos).norm2();     x.real_value += (ai->pos).norm2();
   }   }
 } }
Line 729
Line 708
  
 void colvar::inertia::calc_gradients() void colvar::inertia::calc_gradients()
 { {
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = 2.0 * ai->pos;     ai->grad = 2.0 * ai->pos;
   }   }
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(atoms); 
   } 
 } }
  
  
 void colvar::inertia::apply_force(colvarvalue const &force) void colvar::inertia::apply_force(colvarvalue const &force)
 { {
   if (!atoms.noforce)   if (!atoms->noforce)
     atoms.apply_colvar_force(force.real_value);     atoms->apply_colvar_force(force.real_value);
 } }
  
  
Line 775
Line 749
 void colvar::inertia_z::calc_value() void colvar::inertia_z::calc_value()
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     cvm::real const iprod = ai->pos * axis;     cvm::real const iprod = ai->pos * axis;
     x.real_value += iprod * iprod;     x.real_value += iprod * iprod;
   }   }
Line 784
Line 758
  
 void colvar::inertia_z::calc_gradients() void colvar::inertia_z::calc_gradients()
 { {
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {   for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
     ai->grad = 2.0 * (ai->pos * axis) * axis;     ai->grad = 2.0 * (ai->pos * axis) * axis;
   }   }
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(atoms); 
   } 
 } }
  
  
 void colvar::inertia_z::apply_force(colvarvalue const &force) void colvar::inertia_z::apply_force(colvarvalue const &force)
 { {
   if (!atoms.noforce)   if (!atoms->noforce)
     atoms.apply_colvar_force(force.real_value);     atoms->apply_colvar_force(force.real_value);
 } }
  
  
Line 806
Line 775
 colvar::rmsd::rmsd(std::string const &conf) colvar::rmsd::rmsd(std::string const &conf)
   : cvc(conf)   : cvc(conf)
 { {
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true; 
   function_type = "rmsd";   function_type = "rmsd";
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
   parse_group(conf, "atoms", atoms);   atoms = parse_group(conf, "atoms");
   atom_groups.push_back(&atoms); 
  
   if (atoms.b_dummy) {   if (!atoms || atoms->size() == 0) {
     cvm::error("Error: \"atoms\" cannot be a dummy atom.");     cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD.");
     return;     return;
   }   }
  
   if (atoms.ref_pos_group != NULL && b_Jacobian_derivative) {   bool b_Jacobian_derivative = true;
    if (atoms->fitting_group != NULL && b_Jacobian_derivative) {
     cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "     cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
               "Jacobian derivatives of the RMSD will not be calculated.\n");               "Jacobian derivatives of the RMSD will not be calculated.\n");
     b_Jacobian_derivative = false;     b_Jacobian_derivative = false;
   }   }
    if (b_Jacobian_derivative) provide(f_cvc_Jacobian);
  
   // the following is a simplified version of the corresponding atom group options;   // the following is a simplified version of the corresponding atom group options;
   // we need this because the reference coordinates defined inside the atom group   // we need this because the reference coordinates defined inside the atom group
   // may be used only for fitting, and even more so if ref_pos_group is used   // may be used only for fitting, and even more so if fitting_group is used
   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {   if (get_keyval(conf, "refPositions", ref_pos, ref_pos)) {
     cvm::log("Using reference positions from configuration file to calculate the variable.\n");     cvm::log("Using reference positions from configuration file to calculate the variable.\n");
     if (ref_pos.size() != atoms.size()) {     if (ref_pos.size() != atoms->size()) {
       cvm::error("Error: the number of reference positions provided ("+       cvm::error("Error: the number of reference positions provided ("+
                   cvm::to_str(ref_pos.size())+                   cvm::to_str(ref_pos.size())+
                   ") does not match the number of atoms of group \"atoms\" ("+                   ") does not match the number of atoms of group \"atoms\" ("+
                   cvm::to_str(atoms.size())+").\n");                   cvm::to_str(atoms->size())+").\n");
       return;       return;
     }     }
   }   }
Line 861
Line 830
         }         }
       } else {       } else {
         // if not, rely on existing atom indices for the group         // if not, rely on existing atom indices for the group
         atoms.create_sorted_ids();         atoms->create_sorted_ids();
          ref_pos.resize(atoms->size());
       }       }
  
       ref_pos.resize(atoms.size());       cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms->sorted_ids,
       cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms.sorted_ids, 
                         ref_pos_col, ref_pos_col_value);                         ref_pos_col, ref_pos_col_value);
     }     }
   }   }
  
   if (atoms.b_user_defined_fit) {   if (ref_pos.size() != atoms->size()) {
      cvm::error("Error: found " + cvm::to_str(ref_pos.size()) +
                      " reference positions; expected " + cvm::to_str(atoms->size()));
      return;
    }
  
    if (atoms->b_user_defined_fit) {
     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
   } else {   } else {
     // Default: fit everything     // Default: fit everything
     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
     atoms.b_center = true;     atoms->b_center = true;
     atoms.b_rotate = true;     atoms->b_rotate = true;
     // default case: reference positions for calculating the rmsd are also those used     // default case: reference positions for calculating the rmsd are also those used
     // for fitting     // for fitting
     atoms.ref_pos = ref_pos;     atoms->ref_pos = ref_pos;
     atoms.center_ref_pos();     atoms->center_ref_pos();
  
     cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "     cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
               "will not be computed as they cancel out in the gradients.");               "will not be computed as they cancel out in the gradients.");
     atoms.b_fit_gradients = false;     atoms->b_fit_gradients = false;
   } 
  
   if (atoms.b_rotate) { 
     // TODO: finer-grained control of this would require exposing a 
     // "request_Jacobian_derivative()" method to the colvar, and the same 
     // from the colvar to biases 
     // TODO: this should not be enabled here anyway, as it is not specific of the 
     // component - instead it should be decided in a generic way by the atom group 
  
     // request the calculation of the derivatives of the rotation defined by the atom group     // request the calculation of the derivatives of the rotation defined by the atom group
     atoms.rot.request_group1_gradients(atoms.size());     atoms->rot.request_group1_gradients(atoms->size());
     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
     atoms.rot.request_group2_gradients(atoms.size());     // this is only required for ABF, but we do both groups here for better caching
      atoms->rot.request_group2_gradients(atoms->size());
   }   }
 } }
  
Line 908
Line 876
   // rotational-translational fit is handled by the atom group   // rotational-translational fit is handled by the atom group
  
   x.real_value = 0.0;   x.real_value = 0.0;
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
     x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2();     x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2();
   }   }
   x.real_value /= cvm::real(atoms.size()); // MSD   x.real_value /= cvm::real(atoms->size()); // MSD
   x.real_value = std::sqrt(x.real_value);   x.real_value = std::sqrt(x.real_value);
 } }
  
Line 919
Line 887
 void colvar::rmsd::calc_gradients() void colvar::rmsd::calc_gradients()
 { {
   cvm::real const drmsddx2 = (x.real_value > 0.0) ?   cvm::real const drmsddx2 = (x.real_value > 0.0) ?
     0.5 / (x.real_value * cvm::real(atoms.size())) :     0.5 / (x.real_value * cvm::real(atoms->size())) :
     0.0;     0.0;
  
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
     atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia]));     (*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia]));
   } 
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(atoms); 
   }   }
 } }
  
  
 void colvar::rmsd::apply_force(colvarvalue const &force) void colvar::rmsd::apply_force(colvarvalue const &force)
 { {
   if (!atoms.noforce)   if (!atoms->noforce)
     atoms.apply_colvar_force(force.real_value);     atoms->apply_colvar_force(force.real_value);
 } }
  
  
 void colvar::rmsd::calc_force_invgrads() void colvar::rmsd::calc_force_invgrads()
 { {
   atoms.read_system_forces();   atoms->read_total_forces();
   ft.real_value = 0.0;   ft.real_value = 0.0;
  
   // Note: gradient square norm is 1/N_atoms   // Note: gradient square norm is 1/N_atoms
  
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
     ft.real_value += atoms[ia].grad * atoms[ia].system_force;     ft.real_value += (*atoms)[ia].grad * (*atoms)[ia].total_force;
   }   }
   ft.real_value *= atoms.size();   ft.real_value *= atoms->size();
 } }
  
  
Line 959
Line 922
   // divergence of the rotated coordinates (including only derivatives of the rotation matrix)   // divergence of the rotated coordinates (including only derivatives of the rotation matrix)
   cvm::real divergence = 0.0;   cvm::real divergence = 0.0;
  
   if (atoms.b_rotate) {   if (atoms->b_rotate) {
  
     // gradient of the rotation matrix     // gradient of the rotation matrix
     cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);     cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
     // gradients of products of 2 quaternion components     // gradients of products of 2 quaternion components
     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
     for (size_t ia = 0; ia < atoms.size(); ia++) {     for (size_t ia = 0; ia < atoms->size(); ia++) {
  
       // Gradient of optimal quaternion wrt current Cartesian position       // Gradient of optimal quaternion wrt current Cartesian position
       cvm::vector1d<cvm::rvector> &dq = atoms.rot.dQ0_1[ia];       cvm::vector1d<cvm::rvector> &dq = atoms->rot.dQ0_1[ia];
  
       g11 = 2.0 * (atoms.rot.q)[1]*dq[1];       g11 = 2.0 * (atoms->rot.q)[1]*dq[1];
       g22 = 2.0 * (atoms.rot.q)[2]*dq[2];       g22 = 2.0 * (atoms->rot.q)[2]*dq[2];
       g33 = 2.0 * (atoms.rot.q)[3]*dq[3];       g33 = 2.0 * (atoms->rot.q)[3]*dq[3];
       g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0];       g01 = (atoms->rot.q)[0]*dq[1] + (atoms->rot.q)[1]*dq[0];
       g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0];       g02 = (atoms->rot.q)[0]*dq[2] + (atoms->rot.q)[2]*dq[0];
       g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0];       g03 = (atoms->rot.q)[0]*dq[3] + (atoms->rot.q)[3]*dq[0];
       g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1];       g12 = (atoms->rot.q)[1]*dq[2] + (atoms->rot.q)[2]*dq[1];
       g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1];       g13 = (atoms->rot.q)[1]*dq[3] + (atoms->rot.q)[3]*dq[1];
       g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2];       g23 = (atoms->rot.q)[2]*dq[3] + (atoms->rot.q)[3]*dq[2];
  
       // Gradient of the rotation matrix wrt current Cartesian position       // Gradient of the rotation matrix wrt current Cartesian position
       grad_rot_mat[0][0] = -2.0 * (g22 + g33);       grad_rot_mat[0][0] = -2.0 * (g22 + g33);
Line 1003
Line 966
       }       }
     }     }
   }   }
    jd.real_value = x.real_value > 0.0 ? (3.0 * atoms->size() - 4.0 - divergence) / x.real_value : 0.0;
   jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0; 
 } }
  
  
Line 1013
Line 975
 colvar::eigenvector::eigenvector(std::string const &conf) colvar::eigenvector::eigenvector(std::string const &conf)
   : cvc(conf)   : cvc(conf)
 { {
   b_inverse_gradients = true;   provide(f_cvc_inv_gradient);
   b_Jacobian_derivative = true;   provide(f_cvc_Jacobian);
   function_type = "eigenvector";   function_type = "eigenvector";
   x.type(colvarvalue::type_scalar);   x.type(colvarvalue::type_scalar);
  
   parse_group(conf, "atoms", atoms);   atoms = parse_group(conf, "atoms");
   atom_groups.push_back(&atoms); 
  
  
   {   {
     bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos);     bool const b_inline = get_keyval(conf, "refPositions", ref_pos, ref_pos);
  
     if (b_inline) {     if (b_inline) {
       cvm::log("Using reference positions from input file.\n");       cvm::log("Using reference positions from input file.\n");
       if (ref_pos.size() != atoms.size()) {       if (ref_pos.size() != atoms->size()) {
         cvm::error("Error: reference positions do not "         cvm::error("Error: reference positions do not "
                           "match the number of requested atoms.\n");                           "match the number of requested atoms.\n");
         return;         return;
Line 1054
Line 1014
         }         }
       } else {       } else {
         // if not, use atom indices         // if not, use atom indices
         atoms.create_sorted_ids();         atoms->create_sorted_ids();
       }       }
  
       ref_pos.resize(atoms.size());       ref_pos.resize(atoms->size());
       cvm::load_coords(file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);       cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
     }     }
   }   }
  
    if (ref_pos.size() == 0) {
      cvm::error("Error: reference positions were not provided.\n", INPUT_ERROR);
      return;
    }
  
    if (ref_pos.size() != atoms->size()) {
      cvm::error("Error: reference positions do not "
                 "match the number of requested atoms.\n", INPUT_ERROR);
      return;
    }
  
   // save for later the geometric center of the provided positions (may not be the origin)   // save for later the geometric center of the provided positions (may not be the origin)
   cvm::rvector ref_pos_center(0.0, 0.0, 0.0);   cvm::rvector ref_pos_center(0.0, 0.0, 0.0);
   for (size_t i = 0; i < atoms.size(); i++) {   for (size_t i = 0; i < atoms->size(); i++) {
     ref_pos_center += ref_pos[i];     ref_pos_center += ref_pos[i];
   }   }
   ref_pos_center *= 1.0 / atoms.size();   ref_pos_center *= 1.0 / atoms->size();
  
   if (atoms.b_user_defined_fit) {   if (atoms->b_user_defined_fit) {
     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");     cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
   } else {   } else {
     // default: fit everything     // default: fit everything
     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "     cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");               "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
     atoms.b_center = true;     atoms->b_center = true;
     atoms.b_rotate = true;     atoms->b_rotate = true;
     atoms.ref_pos = ref_pos;     atoms->ref_pos = ref_pos;
     atoms.center_ref_pos();     atoms->center_ref_pos();
      atoms->b_fit_gradients = false; // cancel out if group is fitted on itself
                                      // and cvc is translationally invariant
  
     // request the calculation of the derivatives of the rotation defined by the atom group     // request the calculation of the derivatives of the rotation defined by the atom group
     atoms.rot.request_group1_gradients(atoms.size());     atoms->rot.request_group1_gradients(atoms->size());
     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:     // request derivatives of optimal rotation wrt reference coordinates for Jacobian:
     // this is only required for ABF, but we do both groups here for better caching     // this is only required for ABF, but we do both groups here for better caching
     atoms.rot.request_group2_gradients(atoms.size());     atoms->rot.request_group2_gradients(atoms->size());
   }   }
  
   {   {
Line 1092
Line 1065
     // now load the eigenvector     // now load the eigenvector
     if (b_inline) {     if (b_inline) {
       cvm::log("Using vector components from input file.\n");       cvm::log("Using vector components from input file.\n");
       if (eigenvec.size() != atoms.size()) {       if (eigenvec.size() != atoms->size()) {
         cvm::fatal_error("Error: vector components do not "         cvm::fatal_error("Error: vector components do not "
                           "match the number of requested atoms.\n");                           "match the number of requested atoms->\n");
       }       }
     }     }
  
Line 1117
Line 1090
         }         }
       } else {       } else {
         // if not, use atom indices         // if not, use atom indices
         atoms.create_sorted_ids();         atoms->create_sorted_ids();
       }       }
  
       eigenvec.resize(atoms.size());       eigenvec.resize(atoms->size());
       cvm::load_coords(file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);       cvm::load_coords(file_name.c_str(), eigenvec, atoms->sorted_ids, file_col, file_col_value);
     }     }
   }   }
  
Line 1132
Line 1105
   }   }
  
   cvm::atom_pos eig_center(0.0, 0.0, 0.0);   cvm::atom_pos eig_center(0.0, 0.0, 0.0);
   for (size_t eil = 0; eil < atoms.size(); eil++) {   for (size_t eil = 0; eil < atoms->size(); eil++) {
     eig_center += eigenvec[eil];     eig_center += eigenvec[eil];
   }   }
   eig_center *= 1.0 / atoms.size();   eig_center *= 1.0 / atoms->size();
   cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n");   cvm::log("Geometric center of the provided vector: "+cvm::to_str(eig_center)+"\n");
  
   bool b_difference_vector = false;   bool b_difference_vector = false;
Line 1143
Line 1116
  
   if (b_difference_vector) {   if (b_difference_vector) {
  
     if (atoms.b_center) {     if (atoms->b_center) {
       // both sets should be centered on the origin for fitting       // both sets should be centered on the origin for fitting
       for (size_t i = 0; i < atoms.size(); i++) {       for (size_t i = 0; i < atoms->size(); i++) {
         eigenvec[i] -= eig_center;         eigenvec[i] -= eig_center;
         ref_pos[i]  -= ref_pos_center;         ref_pos[i]  -= ref_pos_center;
       }       }
     }     }
     if (atoms.b_rotate) {     if (atoms->b_rotate) {
       atoms.rot.calc_optimal_rotation(eigenvec, ref_pos);       atoms->rot.calc_optimal_rotation(eigenvec, ref_pos);
       for (size_t i = 0; i < atoms.size(); i++) {       for (size_t i = 0; i < atoms->size(); i++) {
         eigenvec[i] = atoms.rot.rotate(eigenvec[i]);         eigenvec[i] = atoms->rot.rotate(eigenvec[i]);
       }       }
     }     }
     cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");     cvm::log("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
     for (size_t i = 0; i < atoms.size(); i++) {     for (size_t i = 0; i < atoms->size(); i++) {
       eigenvec[i] -= ref_pos[i];       eigenvec[i] -= ref_pos[i];
     }     }
     if (atoms.b_center) {     if (atoms->b_center) {
       // bring back the ref positions to where they were       // bring back the ref positions to where they were
       for (size_t i = 0; i < atoms.size(); i++) {       for (size_t i = 0; i < atoms->size(); i++) {
         ref_pos[i] += ref_pos_center;         ref_pos[i] += ref_pos_center;
       }       }
     }     }
  
   } else {   } else {
     cvm::log("Centering the provided vector to zero.\n");     cvm::log("Centering the provided vector to zero.\n");
     for (size_t i = 0; i < atoms.size(); i++) {     for (size_t i = 0; i < atoms->size(); i++) {
       eigenvec[i] -= eig_center;       eigenvec[i] -= eig_center;
     }     }
   }   }
Line 1178
Line 1151
  
   // for inverse gradients   // for inverse gradients
   eigenvec_invnorm2 = 0.0;   eigenvec_invnorm2 = 0.0;
   for (size_t ein = 0; ein < atoms.size(); ein++) {   for (size_t ein = 0; ein < atoms->size(); ein++) {
     eigenvec_invnorm2 += eigenvec[ein].norm2();     eigenvec_invnorm2 += eigenvec[ein].norm2();
   }   }
   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
  
   if (b_difference_vector) {   if (b_difference_vector) {
     cvm::log("\"differenceVector\" is on: normalizing the vector.\n");     cvm::log("\"differenceVector\" is on: normalizing the vector.\n");
     for (size_t i = 0; i < atoms.size(); i++) {     for (size_t i = 0; i < atoms->size(); i++) {
       eigenvec[i] *= eigenvec_invnorm2;       eigenvec[i] *= eigenvec_invnorm2;
     }     }
   } else {   } else {
Line 1197
Line 1170
 void colvar::eigenvector::calc_value() void colvar::eigenvector::calc_value()
 { {
   x.real_value = 0.0;   x.real_value = 0.0;
   for (size_t i = 0; i < atoms.size(); i++) {   for (size_t i = 0; i < atoms->size(); i++) {
     x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];     x.real_value += ((*atoms)[i].pos - ref_pos[i]) * eigenvec[i];
   }   }
 } }
  
  
 void colvar::eigenvector::calc_gradients() void colvar::eigenvector::calc_gradients()
 { {
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
     atoms[ia].grad = eigenvec[ia];     (*atoms)[ia].grad = eigenvec[ia];
   } 
  
   if (b_debug_gradients) { 
     cvm::log("Debugging gradients:\n"); 
     debug_gradients(atoms); 
   }   }
 } }
  
  
 void colvar::eigenvector::apply_force(colvarvalue const &force) void colvar::eigenvector::apply_force(colvarvalue const &force)
 { {
   if (!atoms.noforce)   if (!atoms->noforce)
     atoms.apply_colvar_force(force.real_value);     atoms->apply_colvar_force(force.real_value);
 } }
  
  
 void colvar::eigenvector::calc_force_invgrads() void colvar::eigenvector::calc_force_invgrads()
 { {
   atoms.read_system_forces();   atoms->read_total_forces();
   ft.real_value = 0.0;   ft.real_value = 0.0;
  
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
     ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *     ft.real_value += eigenvec_invnorm2 * (*atoms)[ia].grad *
       atoms[ia].system_force;       (*atoms)[ia].total_force;
   }   }
 } }
  
Line 1239
Line 1207
 { {
   // gradient of the rotation matrix   // gradient of the rotation matrix
   cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);   cvm::matrix2d<cvm::rvector> grad_rot_mat(3, 3);
   cvm::quaternion &quat0 = atoms.rot.q;   cvm::quaternion &quat0 = atoms->rot.q;
  
   // gradients of products of 2 quaternion components   // gradients of products of 2 quaternion components
   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
  
   cvm::real sum = 0.0;   cvm::real sum = 0.0;
  
   for (size_t ia = 0; ia < atoms.size(); ia++) {   for (size_t ia = 0; ia < atoms->size(); ia++) {
  
     // Gradient of optimal quaternion wrt current Cartesian position     // Gradient of optimal quaternion wrt current Cartesian position
     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t     // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
     // we can just transpose the derivatives of the direct matrix     // we can just transpose the derivatives of the direct matrix
     cvm::vector1d<cvm::rvector> &dq_1 = atoms.rot.dQ0_1[ia];     cvm::vector1d<cvm::rvector> &dq_1 = atoms->rot.dQ0_1[ia];
  
     g11 = 2.0 * quat0[1]*dq_1[1];     g11 = 2.0 * quat0[1]*dq_1[1];
     g22 = 2.0 * quat0[2]*dq_1[2];     g22 = 2.0 * quat0[2]*dq_1[2];
Line 1289
Line 1257
 colvar::cartesian::cartesian(std::string const &conf) colvar::cartesian::cartesian(std::string const &conf)
   : cvc(conf)   : cvc(conf)
 { {
   b_inverse_gradients = false; 
   b_Jacobian_derivative = false; 
   function_type = "cartesian";   function_type = "cartesian";
  
   parse_group(conf, "atoms", atoms);   atoms = parse_group(conf, "atoms");
   atom_groups.push_back(&atoms); 
  
   bool use_x, use_y, use_z;   bool use_x, use_y, use_z;
   get_keyval(conf, "useX", use_x, true);   get_keyval(conf, "useX", use_x, true);
Line 1308
Line 1273
  
   if (axes.size() == 0) {   if (axes.size() == 0) {
     cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");     cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
      return;
   }   }
  
   x.type(colvarvalue::type_vector);   x.type(colvarvalue::type_vector);
   x.vector1d_value.resize(atoms.size() * axes.size());   x.vector1d_value.resize(atoms->size() * axes.size());
 } }
  
  
Line 1319
Line 1285
 { {
   size_t const dim = axes.size();   size_t const dim = axes.size();
   size_t ia, j;   size_t ia, j;
   for (ia = 0; ia < atoms.size(); ia++) {   for (ia = 0; ia < atoms->size(); ia++) {
     for (j = 0; j < dim; j++) {     for (j = 0; j < dim; j++) {
       x.vector1d_value[dim*ia + j] = atoms[ia].pos[axes[j]];       x.vector1d_value[dim*ia + j] = (*atoms)[ia].pos[axes[j]];
     }     }
   }   }
 } }
Line 1339
Line 1305
 { {
   size_t const dim = axes.size();   size_t const dim = axes.size();
   size_t ia, j;   size_t ia, j;
   if (!atoms.noforce) {   if (!atoms->noforce) {
     cvm::rvector f;     cvm::rvector f;
     for (ia = 0; ia < atoms.size(); ia++) {     for (ia = 0; ia < atoms->size(); ia++) {
       for (j = 0; j < dim; j++) {       for (j = 0; j < dim; j++) {
         f[axes[j]] = force.vector1d_value[dim*ia + j];         f[axes[j]] = force.vector1d_value[dim*ia + j];
       }       }
       atoms[ia].apply_force(f);       (*atoms)[ia].apply_force(f);
     }     }
   }   }
 } }

Legend:
Removed in v.1.15 
changed lines
 Added in v.1.16


Made by using version 1.53 of cvs2html