Difference for src/colvarcomp.C from version 1.8 to 1.9

version 1.8version 1.9
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include "colvarmodule.h" #include "colvarmodule.h"
 #include "colvarvalue.h" #include "colvarvalue.h"
Line 10
Line 10
 colvar::cvc::cvc() colvar::cvc::cvc()
   : sup_coeff(1.0),   : sup_coeff(1.0),
     sup_np(1),     sup_np(1),
     b_enabled(true), 
     b_periodic(false),     b_periodic(false),
     b_inverse_gradients(false),     b_try_scalable(true)
     b_Jacobian_derivative(false), {
     b_debug_gradients(false)   init_cvc_requires();
 {} }
  
  
 colvar::cvc::cvc(std::string const &conf) colvar::cvc::cvc(std::string const &conf)
   : sup_coeff(1.0),   : sup_coeff(1.0),
     sup_np(1),     sup_np(1),
     b_enabled(true), 
     b_periodic(false),     b_periodic(false),
     b_inverse_gradients(false),     b_try_scalable(true)
     b_Jacobian_derivative(false), 
     b_debug_gradients(false) 
 { {
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Initializing cvc base object.\n");     cvm::log("Initializing cvc base object.\n");
  
   get_keyval(conf, "name", this->name, std::string(""), parse_silent);   init_cvc_requires();
  
    if (get_keyval(conf, "name", this->name, std::string(""), parse_silent)) {
      // Temporary description until child object is initialized
      description = "cvc " + name;
    } else {
      description = "uninitialized cvc";
    }
  
   get_keyval(conf, "componentCoeff", sup_coeff, 1.0);   get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
   get_keyval(conf, "componentExp", sup_np, 1);   get_keyval(conf, "componentExp", sup_np, 1);
Line 38
Line 41
   get_keyval(conf, "period", period, 0.0);   get_keyval(conf, "period", period, 0.0);
   get_keyval(conf, "wrapAround", wrap_center, 0.0);   get_keyval(conf, "wrapAround", wrap_center, 0.0);
  
   get_keyval(conf, "debugGradients", b_debug_gradients, false, parse_silent);   // All cvcs implement this
    provide(f_cvc_debug_gradient);
    get_keyval_feature((colvarparse *)this, conf, "debugGradients",
                       f_cvc_debug_gradient, false, parse_silent);
  
    // Attempt scalable calculations when in parallel? (By default yes, if available)
    get_keyval(conf, "scalable", b_try_scalable, true);
  
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Done initializing cvc base object.\n");     cvm::log("Done initializing cvc base object.\n");
 } }
  
  
 void colvar::cvc::parse_group(std::string const &conf, int colvar::cvc::init_total_force_params(std::string const &conf)
  {
    if (get_keyval_feature(this, conf, "oneSiteSystemForce",
                           f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
      cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
               "please use \"oneSiteTotalForce\" instead.\n");
    }
    if (get_keyval_feature(this, conf, "oneSiteTotalForce",
                           f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
      cvm::log("Computing total force on group 1 only");
    }
    return COLVARS_OK;
  }
  
  
  cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
                                char const *group_key,                                char const *group_key,
                                cvm::atom_group &group, 
                                bool optional)                                bool optional)
 { {
    cvm::atom_group *group = NULL;
  
   if (key_lookup(conf, group_key)) {   if (key_lookup(conf, group_key)) {
     if (group.parse(conf, group_key) != COLVARS_OK) {     group = new cvm::atom_group;
      group->key = group_key;
  
      if (b_try_scalable) {
        // TODO rewrite this logic in terms of dependencies
        if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
          enable(f_cvc_scalable_com);
          enable(f_cvc_scalable);
          group->enable(f_ag_scalable_com);
          group->enable(f_ag_scalable);
        }
  
        // TODO check for other types of parallelism here
  
        if (is_enabled(f_cvc_scalable)) {
          cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
        }
      }
  
      if (group->parse(conf) == COLVARS_OK) {
        atom_groups.push_back(group);
      } else {
       cvm::error("Error parsing definition for atom group \""+       cvm::error("Error parsing definition for atom group \""+
                          std::string(group_key)+"\".\n");                          std::string(group_key)+"\".\n");
       return; 
     }     }
   } else {   } else {
     if (! optional) {     if (! optional) {
       cvm::error("Error: definition for atom group \""+       cvm::error("Error: definition for atom group \""+
                       std::string(group_key)+"\" not found.\n");                       std::string(group_key)+"\" not found.\n");
       return; 
     }     }
   }   }
    return group;
  }
  
  
  int colvar::cvc::setup()
  {
    size_t i;
    description = "cvc " + name;
  
    for (i = 0; i < atom_groups.size(); i++) {
      add_child((colvardeps *) atom_groups[i]);
    }
  
    return COLVARS_OK;
 } }
  
  
 colvar::cvc::~cvc() colvar::cvc::~cvc()
 {} {
    remove_all_children();
    for (size_t i = 0; i < atom_groups.size(); i++) {
      if (atom_groups[i] != NULL) delete atom_groups[i];
    }
  }
  
  
  void colvar::cvc::read_data()
  {
    size_t ig;
    for (ig = 0; ig < atom_groups.size(); ig++) {
      cvm::atom_group &atoms = *(atom_groups[ig]);
      atoms.reset_atoms_data();
      atoms.read_positions();
      atoms.calc_required_properties();
      // each atom group will take care of its own fitting_group, if defined
    }
  
  ////  Don't try to get atom velocities, as no back-end currently implements it
  //   if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
  //     for (i = 0; i < cvcs.size(); i++) {
  //       for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
  //         cvcs[i]->atom_groups[ig]->read_velocities();
  //       }
  //     }
  //   }
  }
  
  
 void colvar::cvc::calc_force_invgrads() void colvar::cvc::calc_force_invgrads()
Line 84
Line 169
 } }
  
  
 void colvar::cvc::debug_gradients(cvm::atom_group &group) void colvar::cvc::debug_gradients(cvm::atom_group *group)
 { {
   // this function should work for any scalar variable:   // this function should work for any scalar variable:
   // the only difference will be the name of the atom group (here, "group")   // the only difference will be the name of the atom group (here, "group")
    // NOTE: this assumes that groups for this cvc are non-overlapping,
    // since atom coordinates are modified only within the current group
  
   if (group.b_dummy) return;   if (group->b_dummy) return;
  
   cvm::rotation const rot_0 = group.rot;   cvm::rotation const rot_0 = group->rot;
   cvm::rotation const rot_inv = group.rot.inverse();   cvm::rotation const rot_inv = group->rot.inverse();
  
   cvm::real x_0 = x.real_value;   cvm::real x_0 = x.real_value;
   if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0];   if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0];
  
   // cvm::log("gradients     = "+cvm::to_str (gradients)+"\n");   // cvm::log("gradients     = "+cvm::to_str (gradients)+"\n");
  
   // it only makes sense to debug the fit gradients   cvm::atom_group *group_for_fit = group->fitting_group ? group->fitting_group : group;
   // when the fitting group is the same as this group   cvm::atom_pos fit_gradient_sum, gradient_sum;
   if (group.b_rotate || group.b_center) 
     if (group.b_fit_gradients && (group.ref_pos_group == NULL)) {   // print the values of the fit gradients
       group.calc_fit_gradients();   if (group->b_rotate || group->b_center) {
       if (group.b_rotate) {     if (group->b_fit_gradients) {
         // fit_gradients are in the original frame, we should print them in the rotated frame       size_t j;
         for (size_t j = 0; j < group.fit_gradients.size(); j++) { 
           group.fit_gradients[j] = rot_0.rotate(group.fit_gradients[j]);       // fit_gradients are in the simulation frame: we should print them in the rotated frame
         }       cvm::log("Fit gradients:\n");
       }       for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
       cvm::log("fit_gradients = "+cvm::to_str(group.fit_gradients)+"\n");         cvm::log((group->fitting_group ? std::string("refPosGroup") : group->key) +
       if (group.b_rotate) {                  "[" + cvm::to_str(j) + "] = " +
         for (size_t j = 0; j < group.fit_gradients.size(); j++) {                  (group->b_rotate ?
           group.fit_gradients[j] = rot_inv.rotate(group.fit_gradients[j]);                   cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) :
                    cvm::to_str(group_for_fit->fit_gradients[j])));
         }         }
       }       }
     }     }
  
   for (size_t ia = 0; ia < group.size(); ia++) {   // debug the gradients
    for (size_t ia = 0; ia < group->size(); ia++) {
  
     // tests are best conducted in the unrotated (simulation) frame     // tests are best conducted in the unrotated (simulation) frame
     cvm::rvector const atom_grad = group.b_rotate ?     cvm::rvector const atom_grad = (group->b_rotate ?
       rot_inv.rotate(group[ia].grad) :                                     rot_inv.rotate((*group)[ia].grad) :
       group[ia].grad;                                     (*group)[ia].grad);
      gradient_sum += atom_grad;
  
     for (size_t id = 0; id < 3; id++) {     for (size_t id = 0; id < 3; id++) {
       // (re)read original positions       // (re)read original positions
       group.read_positions();       group->read_positions();
       // change one coordinate       // change one coordinate
       group[ia].pos[id] += cvm::debug_gradients_step_size;       (*group)[ia].pos[id] += cvm::debug_gradients_step_size;
       // (re)do the fit (if defined)       group->calc_required_properties();
       if (group.b_center || group.b_rotate) { 
         group.calc_apply_roto_translation(); 
       } 
       calc_value();       calc_value();
       cvm::real x_1 = x.real_value;       cvm::real x_1 = x.real_value;
       if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];       if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
       cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n");       cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n");
       cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0,       cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0,
                              21, 14)+"\n");                              21, 14)+"\n");
       //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ?       cvm::real const dx_pred = (group->fit_gradients.size()) ?
       cvm::real const dx_pred = (group.fit_gradients.size()) ?         (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
         (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : 
         (cvm::debug_gradients_step_size * atom_grad[id]);         (cvm::debug_gradients_step_size * atom_grad[id]);
       cvm::log("dx(interp) = "+cvm::to_str(dx_pred,       cvm::log("dx(interp) = "+cvm::to_str(dx_pred,
                              21, 14)+"\n");                              21, 14)+"\n");
Line 152
Line 238
     }     }
   }   }
  
 /*   if ((group->b_fit_gradients) && (group->fitting_group != NULL)) {
  * The code below is WIP     cvm::atom_group *ref_group = group->fitting_group;
  */     group->read_positions();
 //   if (group.ref_pos_group != NULL) {     group->calc_required_properties();
 //     cvm::atom_group &ref = *group.ref_pos_group; 
 //     group.calc_fit_gradients();     for (size_t ia = 0; ia < ref_group->size(); ia++) {
 // 
 //     for (size_t ia = 0; ia < ref.size(); ia++) {       // fit gradients are in the unrotated (simulation) frame
 //       cvm::rvector const atom_grad = ref_group->fit_gradients[ia];
 //       for (size_t id = 0; id < 3; id++) {       fit_gradient_sum += atom_grad;
 //         // (re)read original positions 
 //         group.read_positions();       for (size_t id = 0; id < 3; id++) {
 //         ref.read_positions();         // (re)read original positions
 //         // change one coordinate         group->read_positions();
 //         ref[ia].pos[id] += cvm::debug_gradients_step_size;         ref_group->read_positions();
 //         group.calc_apply_roto_translation();         // change one coordinate
 //         calc_value();         (*ref_group)[ia].pos[id] += cvm::debug_gradients_step_size;
 //         cvm::real const x_1 = x.real_value;         group->calc_required_properties();
 //         cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");         calc_value();
 //         cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0, 
 //                                21, 14)+"\n");         cvm::real const x_1 = x.real_value;
 //         //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ?         cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");
 //         // cvm::real const dx_pred = (group.fit_gradients.size()) ?         cvm::log("dx(actual) = "+cvm::to_str (x_1 - x_0,
 //         //  (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) :                                21, 14)+"\n");
 //         //  (cvm::debug_gradients_step_size * atom_grad[id]); 
 //         cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][id];         cvm::real const dx_pred = cvm::debug_gradients_step_size * atom_grad[id];
 //         cvm::log("dx(interp) = "+cvm::to_str (dx_pred, 
 //                                21, 14)+"\n");         cvm::log("dx(interp) = "+cvm::to_str (dx_pred,
 //         cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+                                21, 14)+"\n");
 //                   cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /         cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
 //                                std::fabs (x_1 - x_0),                   cvm::to_str(std::fabs (x_1 - x_0 - dx_pred) /
 //                                12, 5)+                                std::fabs (x_1 - x_0),
 //                   ".\n");                                12, 5)+
 //       }                   ".\n");
 //     }       }
 //   }     }
    }
  
    cvm::log("Gradient sum: " +  cvm::to_str(gradient_sum) +
          "  Fit gradient sum: " + cvm::to_str(fit_gradient_sum) +
          "  Total " + cvm::to_str(gradient_sum + fit_gradient_sum));
  
   return;   return;
 } }
  
  
  // Static members
  
  std::vector<colvardeps::feature *> colvar::cvc::cvc_features;

Legend:
Removed in v.1.8 
changed lines
 Added in v.1.9


Made by using version 1.53 of cvs2html