diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 74e52f507..40f3c236c 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -1,1615 +1,1672 @@ // -*- c++ -*- #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" #include "colvar.h" #include "colvarcomp.h" #include // XX TODO make the acf code handle forces as well as values and velocities colvar::colvar (std::string const &conf) { cvm::log ("Initializing a new collective variable.\n"); get_keyval (conf, "name", this->name, (std::string ("colvar")+cvm::to_str (cvm::colvars.size()+1))); for (std::vector::iterator cvi = cvm::colvars.begin(); cvi < cvm::colvars.end(); cvi++) { if ((*cvi)->name == this->name) cvm::fatal_error ("Error: this colvar cannot have the same name, \""+this->name+ "\", as another colvar.\n"); } // all tasks disabled by default for (size_t i = 0; i < task_ntot; i++) { tasks[i] = false; } kinetic_energy = 0.0; potential_energy = 0.0; // read the configuration and set up corresponding instances, for // each type of component implemented #define initialize_components(def_desc,def_config_key,def_class_name) \ { \ size_t def_count = 0; \ std::string def_conf = ""; \ size_t pos = 0; \ while ( this->key_lookup (conf, \ def_config_key, \ def_conf, \ pos) ) { \ if (!def_conf.size()) continue; \ cvm::log ("Initializing " \ "a new \""+std::string (def_config_key)+"\" component"+ \ (cvm::debug() ? ", with configuration:\n"+def_conf \ : ".\n")); \ cvm::increase_depth(); \ cvc *cvcp = new colvar::def_class_name (def_conf); \ if (cvcp != NULL) { \ cvcs.push_back (cvcp); \ cvcp->check_keywords (def_conf, def_config_key); \ cvm::decrease_depth(); \ } else { \ cvm::fatal_error ("Error: in allocating component \"" \ def_config_key"\".\n"); \ } \ if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) { \ if ( (cvcp->function_type != std::string ("distance_z")) && \ (cvcp->function_type != std::string ("dihedral")) && \ (cvcp->function_type != std::string ("spin_angle")) ) { \ cvm::fatal_error ("Error: invalid use of period and/or " \ "wrapAround in a \""+ \ std::string (def_config_key)+ \ "\" component.\n"+ \ "Period: "+cvm::to_str(cvcp->period) + \ " wrapAround: "+cvm::to_str(cvcp->wrap_center));\ } \ } \ if ( ! cvcs.back()->name.size()) \ cvcs.back()->name = std::string (def_config_key)+ \ (cvm::to_str (++def_count)); \ if (cvm::debug()) \ cvm::log ("Done initializing a \""+ \ std::string (def_config_key)+ \ "\" component"+ \ (cvm::debug() ? \ ", named \""+cvcs.back()->name+"\"" \ : "")+".\n"); \ def_conf = ""; \ if (cvm::debug()) \ cvm::log ("Parsed "+cvm::to_str (cvcs.size())+ \ " components at this time.\n"); \ } \ } initialize_components ("distance", "distance", distance); initialize_components ("distance vector", "distanceVec", distance_vec); initialize_components ("distance vector " "direction", "distanceDir", distance_dir); initialize_components ("distance projection " "on an axis", "distanceZ", distance_z); initialize_components ("distance projection " "on a plane", "distanceXY", distance_xy); initialize_components ("average distance weighted by inverse power", "distanceInv", distance_inv); initialize_components ("coordination " "number", "coordNum", coordnum); initialize_components ("self-coordination " "number", "selfCoordNum", selfcoordnum); initialize_components ("angle", "angle", angle); initialize_components ("dihedral", "dihedral", dihedral); initialize_components ("hydrogen bond", "hBond", h_bond); // initialize_components ("alpha helix", "alphaDihedrals", alpha_dihedrals); initialize_components ("alpha helix", "alpha", alpha_angles); initialize_components ("dihedral principal " "component", "dihedralPC", dihedPC); initialize_components ("orientation", "orientation", orientation); initialize_components ("orientation " "angle", "orientationAngle",orientation_angle); initialize_components ("tilt", "tilt", tilt); initialize_components ("spin angle", "spinAngle", spin_angle); initialize_components ("RMSD", "rmsd", rmsd); // initialize_components ("logarithm of MSD", "logmsd", logmsd); initialize_components ("radius of " "gyration", "gyration", gyration); initialize_components ("moment of " "inertia", "inertia", inertia); initialize_components ("moment of inertia around an axis", "inertiaZ", inertia_z); initialize_components ("eigenvector", "eigenvector", eigenvector); if (!cvcs.size()) cvm::fatal_error ("Error: no valid components were provided " "for this collective variable.\n"); cvm::log ("All components initialized.\n"); // this is set false if any of the components has an exponent // different from 1 in the polynomial b_linear = true; // these will be set to false if any of the cvcs has them false b_inverse_gradients = true; b_Jacobian_force = true; // Decide whether the colvar is periodic // Used to wrap extended DOF if extendedLagrangian is on if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1 && (cvcs[0])->sup_coeff == 1.0 ) { this->b_periodic = true; this->period = (cvcs[0])->period; // TODO write explicit wrap() function for colvars to allow for // sup_coeff different from 1 // this->period = (cvcs[0])->period * (cvcs[0])->sup_coeff; } else { this->b_periodic = false; this->period = 0.0; } // check the available features of each cvc for (size_t i = 0; i < cvcs.size(); i++) { if ((cvcs[i])->sup_np != 1) { if (cvm::debug() && b_linear) cvm::log ("Warning: You are using a non-linear polynomial " "combination to define this collective variable, " "some biasing methods may be unavailable.\n"); b_linear = false; if ((cvcs[i])->sup_np < 0) { cvm::log ("Warning: you chose a negative exponent in the combination; " "if you apply forces, the simulation may become unstable " "when the component \""+ (cvcs[i])->function_type+"\" approaches zero.\n"); } } if ((cvcs[i])->b_periodic && !b_periodic) { cvm::log ("Warning: although this component is periodic, the colvar will " "not be treated as periodic, either because the exponent is not " "1, or because multiple components are present. Make sure that " "you know what you are doing!"); } if (! (cvcs[i])->b_inverse_gradients) b_inverse_gradients = false; if (! (cvcs[i])->b_Jacobian_derivative) b_Jacobian_force = false; for (size_t j = i; j < cvcs.size(); j++) { if ( (cvcs[i])->type() != (cvcs[j])->type() ) { cvm::fatal_error ("ERROR: you are definining this collective variable " "by using components of different types, \""+ colvarvalue::type_desc[(cvcs[i])->type()]+ "\" and \""+ colvarvalue::type_desc[(cvcs[j])->type()]+ "\". " "You must use the same type in order to " " sum them together.\n"); } } } { colvarvalue::Type const value_type = (cvcs[0])->type(); if (cvm::debug()) cvm::log ("This collective variable is a "+ colvarvalue::type_desc[value_type]+", corresponding to "+ cvm::to_str (colvarvalue::dof_num[value_type])+ " internal degrees of freedom.\n"); x.type (value_type); x_reported.type (value_type); } get_keyval (conf, "width", width, 1.0); if (width <= 0.0) cvm::fatal_error ("Error: \"width\" must be positive.\n"); lower_boundary.type (this->type()); lower_wall.type (this->type()); upper_boundary.type (this->type()); upper_wall.type (this->type()); if (this->type() == colvarvalue::type_scalar) { if (get_keyval (conf, "lowerBoundary", lower_boundary, lower_boundary)) { enable (task_lower_boundary); } get_keyval (conf, "lowerWallConstant", lower_wall_k, 0.0); if (lower_wall_k > 0.0) { get_keyval (conf, "lowerWall", lower_wall, lower_boundary); enable (task_lower_wall); } if (get_keyval (conf, "upperBoundary", upper_boundary, upper_boundary)) { enable (task_upper_boundary); } get_keyval (conf, "upperWallConstant", upper_wall_k, 0.0); if (upper_wall_k > 0.0) { get_keyval (conf, "upperWall", upper_wall, upper_boundary); enable (task_upper_wall); } } if (tasks[task_lower_boundary]) { get_keyval (conf, "hardLowerBoundary", hard_lower_boundary, false); } if (tasks[task_upper_boundary]) { get_keyval (conf, "hardUpperBoundary", hard_upper_boundary, false); } // consistency checks for boundaries and walls if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) { if (lower_boundary >= upper_boundary) { cvm::fatal_error ("Error: the upper boundary, "+ cvm::to_str (upper_boundary)+ ", is not higher than the lower boundary, "+ cvm::to_str (lower_boundary)+".\n"); } } if (tasks[task_lower_wall] && tasks[task_upper_wall]) { if (lower_wall >= upper_wall) { cvm::fatal_error ("Error: the upper wall, "+ cvm::to_str (upper_wall)+ ", is not higher than the lower wall, "+ cvm::to_str (lower_wall)+".\n"); } if (dist2 (lower_wall, upper_wall) < 1.0E-12) { cvm::log ("Lower wall and upper wall are equal " "in the periodic domain of the colvar: disabling walls.\n"); disable (task_lower_wall); disable (task_upper_wall); } } get_keyval (conf, "expandBoundaries", expand_boundaries, false); if (expand_boundaries && periodic_boundaries()) { cvm::fatal_error ("Error: trying to expand boundaries that already " "cover a whole period of a periodic colvar.\n"); } + if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) { + cvm::fatal_error ("Error: inconsistent configuration " + "(trying to expand boundaries with both " + "hardLowerBoundary and hardUpperBoundary enabled).\n"); + } { bool b_extended_lagrangian; get_keyval (conf, "extendedLagrangian", b_extended_lagrangian, false); if (b_extended_lagrangian) { cvm::real temp, tolerance, period; cvm::log ("Enabling the extended Lagrangian term for colvar \""+ this->name+"\".\n"); enable (task_extended_lagrangian); xr.type (this->type()); vr.type (this->type()); fr.type (this->type()); const bool found = get_keyval (conf, "extendedTemp", temp, cvm::temperature()); if (temp <= 0.0) { if (found) cvm::fatal_error ("Error: \"extendedTemp\" must be positive.\n"); else cvm::fatal_error ("Error: a positive temperature must be provided, either " "by enabling a thermostat, or through \"extendedTemp\".\n"); } get_keyval (conf, "extendedFluctuation", tolerance, width); if (tolerance <= 0.0) cvm::fatal_error ("Error: \"extendedFluctuation\" must be positive.\n"); ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); cvm::log ("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2"); get_keyval (conf, "extendedTimeConstant", period, 200.0); if (period <= 0.0) cvm::fatal_error ("Error: \"extendedTimeConstant\" must be positive.\n"); ext_mass = (cvm::boltzmann() * temp * period * period) / (4.0 * PI * PI * tolerance * tolerance); cvm::log ("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)"); { bool b_output_energy; get_keyval (conf, "outputEnergy", b_output_energy, false); if (b_output_energy) { enable (task_output_energy); } } get_keyval (conf, "extendedLangevinDamping", ext_gamma, 1.0); if (ext_gamma < 0.0) cvm::fatal_error ("Error: \"extendedLangevinDamping\" may not be negative.\n"); if (ext_gamma != 0.0) { enable (task_langevin); ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1 ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt()); } } } { bool b_output_value; get_keyval (conf, "outputValue", b_output_value, true); if (b_output_value) { enable (task_output_value); } } { bool b_output_velocity; get_keyval (conf, "outputVelocity", b_output_velocity, false); if (b_output_velocity) { enable (task_output_velocity); } } { bool b_output_system_force; get_keyval (conf, "outputSystemForce", b_output_system_force, false); if (b_output_system_force) { enable (task_output_system_force); } } { bool b_output_applied_force; get_keyval (conf, "outputAppliedForce", b_output_applied_force, false); if (b_output_applied_force) { enable (task_output_applied_force); } } if (cvm::b_analysis) parse_analysis (conf); if (cvm::debug()) cvm::log ("Done initializing collective variable \""+this->name+"\".\n"); } void colvar::build_atom_list (void) { // If atomic gradients are requested, build full list of atom ids from all cvcs std::list temp_id_list; for (size_t i = 0; i < cvcs.size(); i++) { for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { temp_id_list.push_back (cvcs[i]->atom_groups[j]->at(k).id); } } } temp_id_list.sort(); temp_id_list.unique(); atom_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); temp_id_list.clear(); atomic_gradients.resize (atom_ids.size()); if (atom_ids.size()) { if (cvm::debug()) cvm::log ("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n"); } else { cvm::log ("Warning: colvar components communicated no atom IDs.\n"); } } void colvar::parse_analysis (std::string const &conf) { // if (cvm::debug()) // cvm::log ("Parsing analysis flags for collective variable \""+ // this->name+"\".\n"); runave_length = 0; bool b_runave = false; if (get_keyval (conf, "runAve", b_runave) && b_runave) { enable (task_runave); get_keyval (conf, "runAveLength", runave_length, 1000); get_keyval (conf, "runAveStride", runave_stride, 1); if ((cvm::restart_out_freq % runave_stride) != 0) cvm::fatal_error ("Error: runAveStride must be commensurate with the restart frequency.\n"); std::string runave_outfile; get_keyval (conf, "runAveOutputFile", runave_outfile, std::string (cvm::output_prefix+"."+ this->name+".runave.traj")); size_t const this_cv_width = x.output_width (cvm::cv_width); runave_os.open (runave_outfile.c_str()); runave_os << "# " << cvm::wrap_string ("step", cvm::it_width-2) << " " << cvm::wrap_string ("running average", this_cv_width) << " " << cvm::wrap_string ("running stddev", this_cv_width) << "\n"; } acf_length = 0; bool b_acf = false; if (get_keyval (conf, "corrFunc", b_acf) && b_acf) { enable (task_corrfunc); std::string acf_colvar_name; get_keyval (conf, "corrFuncWithColvar", acf_colvar_name, this->name); if (acf_colvar_name == this->name) { cvm::log ("Calculating auto-correlation function.\n"); } else { cvm::log ("Calculating correlation function with \""+ this->name+"\".\n"); } std::string acf_type_str; get_keyval (conf, "corrFuncType", acf_type_str, to_lower_cppstr (std::string ("velocity"))); if (acf_type_str == to_lower_cppstr (std::string ("coordinate"))) { acf_type = acf_coor; } else if (acf_type_str == to_lower_cppstr (std::string ("velocity"))) { acf_type = acf_vel; enable (task_fdiff_velocity); if (acf_colvar_name.size()) (cvm::colvar_p (acf_colvar_name))->enable (task_fdiff_velocity); } else if (acf_type_str == to_lower_cppstr (std::string ("coordinate_p2"))) { acf_type = acf_p2coor; } else { cvm::fatal_error ("Unknown type of correlation function, \""+ acf_type_str+"\".\n"); } get_keyval (conf, "corrFuncOffset", acf_offset, 0); get_keyval (conf, "corrFuncLength", acf_length, 1000); get_keyval (conf, "corrFuncStride", acf_stride, 1); if ((cvm::restart_out_freq % acf_stride) != 0) cvm::fatal_error ("Error: corrFuncStride must be commensurate with the restart frequency.\n"); get_keyval (conf, "corrFuncNormalize", acf_normalize, true); get_keyval (conf, "corrFuncOutputFile", acf_outfile, std::string (cvm::output_prefix+"."+this->name+ ".corrfunc.dat")); } } void colvar::enable (colvar::task const &t) { switch (t) { case task_output_system_force: enable (task_system_force); break; case task_report_Jacobian_force: enable (task_Jacobian_force); enable (task_system_force); if (cvm::debug()) cvm::log ("Adding the Jacobian force to the system force, " "rather than applying its opposite silently.\n"); break; case task_Jacobian_force: // checks below do not apply to extended-system colvars if ( !tasks[task_extended_lagrangian] ) { enable (task_gradients); if (!b_Jacobian_force) cvm::fatal_error ("Error: colvar \""+this->name+ "\" does not have Jacobian forces implemented.\n"); if (!b_linear) cvm::fatal_error ("Error: colvar \""+this->name+ "\" must be defined as a linear combination " "to calculate the Jacobian force.\n"); if (cvm::debug()) cvm::log ("Enabling calculation of the Jacobian force " "on this colvar.\n"); } fj.type (this->type()); break; case task_system_force: if (!tasks[task_extended_lagrangian]) { if (!b_inverse_gradients) { cvm::fatal_error ("Error: one or more of the components of " "colvar \""+this->name+ "\" does not support system force calculation.\n"); } cvm::request_system_force(); } ft.type (this->type()); ft_reported.type (this->type()); break; case task_output_applied_force: case task_lower_wall: case task_upper_wall: // all of the above require gradients enable (task_gradients); break; case task_fdiff_velocity: x_old.type (this->type()); v_fdiff.type (this->type()); v_reported.type (this->type()); break; case task_output_velocity: enable (task_fdiff_velocity); break; case task_grid: if (this->type() != colvarvalue::type_scalar) cvm::fatal_error ("Cannot calculate a grid for collective variable, \""+ this->name+"\", because its value is not a scalar number.\n"); break; case task_extended_lagrangian: enable (task_gradients); v_reported.type (this->type()); break; case task_lower_boundary: case task_upper_boundary: if (this->type() != colvarvalue::type_scalar) { cvm::fatal_error ("Error: this colvar is not a scalar value " "and cannot produce a grid.\n"); } break; case task_output_value: case task_runave: case task_corrfunc: case task_ntot: break; case task_gradients: f.type (this->type()); fb.type (this->type()); break; case task_collect_gradients: if (this->type() != colvarvalue::type_scalar) cvm::fatal_error ("Collecting atomic gradients for non-scalar collective variable \""+ this->name+"\" is not yet implemented.\n"); enable (task_gradients); if (atom_ids.size() == 0) { build_atom_list(); } break; } tasks[t] = true; } void colvar::disable (colvar::task const &t) { // check dependencies switch (t) { case task_gradients: disable (task_upper_wall); disable (task_lower_wall); disable (task_output_applied_force); disable (task_system_force); disable (task_Jacobian_force); break; case task_system_force: disable (task_output_system_force); break; case task_Jacobian_force: disable (task_report_Jacobian_force); break; case task_fdiff_velocity: disable (task_output_velocity); break; case task_lower_boundary: case task_upper_boundary: disable (task_grid); break; case task_extended_lagrangian: case task_report_Jacobian_force: case task_output_value: case task_output_velocity: case task_output_applied_force: case task_output_system_force: case task_runave: case task_corrfunc: case task_grid: case task_lower_wall: case task_upper_wall: case task_ntot: break; } tasks[t] = false; } colvar::~colvar() { for (size_t i = 0; i < cvcs.size(); i++) { delete cvcs[i]; } } // ******************** CALC FUNCTIONS ******************** void colvar::calc() { if (cvm::debug()) cvm::log ("Calculating colvar \""+this->name+"\".\n"); + // prepare atom groups for calculation + for (size_t i = 0; i < cvcs.size(); i++) { + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); + atoms.reset_atoms_data(); + atoms.read_positions(); + if (atoms.b_center || atoms.b_rotate) { + atoms.calc_apply_roto_translation(); + } + // each atom group will take care of its own ref_pos_group, if defined + } + } + if (tasks[task_output_velocity]) { + for (size_t i = 0; i < cvcs.size(); i++) { + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + cvcs[i]->atom_groups[ig]->read_velocities(); + } + } + } + if (tasks[task_system_force]) { + for (size_t i = 0; i < cvcs.size(); i++) { + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + cvcs[i]->atom_groups[ig]->read_system_forces(); + } + } + } + // calculate the value of the colvar x.reset(); if (x.type() == colvarvalue::type_scalar) { + // polynomial combination allowed - // scalar variable, polynomial combination allowed for (size_t i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->calc_value(); cvm::decrease_depth(); if (cvm::debug()) cvm::log ("Colvar component no. "+cvm::to_str (i+1)+ " within colvar \""+this->name+"\" has value "+ cvm::to_str ((cvcs[i])->value(), cvm::cv_width, cvm::cv_prec)+".\n"); x += (cvcs[i])->sup_coeff * ( ((cvcs[i])->sup_np != 1) ? std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : (cvcs[i])->value().real_value ); } } else { + // only linear combination allowed for (size_t i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->calc_value(); cvm::decrease_depth(); if (cvm::debug()) cvm::log ("Colvar component no. "+cvm::to_str (i+1)+ " within colvar \""+this->name+"\" has value "+ cvm::to_str ((cvcs[i])->value(), cvm::cv_width, cvm::cv_prec)+".\n"); x += (cvcs[i])->sup_coeff * (cvcs[i])->value(); } } if (cvm::debug()) cvm::log ("Colvar \""+this->name+"\" has value "+ cvm::to_str (x, cvm::cv_width, cvm::cv_prec)+".\n"); if (tasks[task_gradients]) { - // calculate the gradients + + if (cvm::debug()) + cvm::log ("Calculating gradients of colvar \""+this->name+"\".\n"); + for (size_t i = 0; i < cvcs.size(); i++) { + // calculate the gradients of each component cvm::increase_depth(); + (cvcs[i])->calc_gradients(); + + // if requested, propagate (via chain rule) the gradients above + // to the atoms used to define the roto-translation + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + if (cvcs[i]->atom_groups[ig]->b_fit_gradients) + cvcs[i]->atom_groups[ig]->calc_fit_gradients(); + } + cvm::decrease_depth(); } + + if (cvm::debug()) + cvm::log ("Done calculating gradients of colvar \""+this->name+"\".\n"); + if (tasks[task_collect_gradients]) { // Collect the atomic gradients inside colvar object for (int a = 0; a < atomic_gradients.size(); a++) { atomic_gradients[a].reset(); } for (size_t i = 0; i < cvcs.size(); i++) { // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) * std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { // If necessary, apply inverse rotation to get atomic // gradient in the laboratory frame if (cvcs[i]->atom_groups[j]->b_rotate) { cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse(); for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { int a = std::lower_bound (atom_ids.begin(), atom_ids.end(), cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin(); atomic_gradients[a] += coeff * rot_inv.rotate (cvcs[i]->atom_groups[j]->at(k).grad); } } else { for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { int a = std::lower_bound (atom_ids.begin(), atom_ids.end(), cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin(); atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad; } } } } } } if (tasks[task_system_force]) { + + if (cvm::debug()) + cvm::log ("Calculating system force of colvar \""+this->name+"\".\n"); + ft.reset(); if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { // get from the cvcs the system forces from the PREVIOUS step for (size_t i = 0; i < cvcs.size(); i++) { (cvcs[i])->calc_force_invgrads(); // linear combination is assumed cvm::increase_depth(); ft += (cvcs[i])->system_force() / ((cvcs[i])->sup_coeff * cvm::real (cvcs.size())); cvm::decrease_depth(); } } if (tasks[task_report_Jacobian_force]) { // add the Jacobian force to the system force, and don't apply any silent // correction internally: biases such as colvarbias_abf will handle it ft += fj; } + + if (cvm::debug()) + cvm::log ("Done calculating system force of colvar \""+this->name+"\".\n"); } if (tasks[task_fdiff_velocity]) { // calculate the velocity by finite differences if (cvm::step_relative() == 0) x_old = x; else { v_fdiff = fdiff_velocity (x_old, x); v_reported = v_fdiff; } } if (tasks[task_extended_lagrangian]) { // initialize the restraint center in the first step to the value // just calculated from the cvcs // TODO: put it in the restart information if (cvm::step_relative() == 0) { xr = x; vr = 0.0; // (already 0; added for clarity) } // report the restraint center as "value" x_reported = xr; v_reported = vr; // the "system force" with the extended Lagrangian is just the // harmonic term acting on the extended coordinate // Note: this is the force for current timestep ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad (xr, x); } else { x_reported = x; ft_reported = ft; } if (cvm::debug()) cvm::log ("Done calculating colvar \""+this->name+"\".\n"); } cvm::real colvar::update() { if (cvm::debug()) cvm::log ("Updating colvar \""+this->name+"\".\n"); // set to zero the applied force f.reset(); // add the biases' force, which at this point should already have // been summed over each bias using this colvar f += fb; if (tasks[task_lower_wall] || tasks[task_upper_wall]) { // wall force colvarvalue fw (this->type()); // if the two walls are applied concurrently, decide which is the // closer one (on a periodic colvar, both walls may be applicable // at the same time) if ( (!tasks[task_upper_wall]) || (this->dist2 (x, lower_wall) < this->dist2 (x, upper_wall)) ) { cvm::real const grad = this->dist2_lgrad (x, lower_wall); if (grad < 0.0) { fw = -0.5 * lower_wall_k * grad; if (cvm::debug()) cvm::log ("Applying a lower wall force ("+ cvm::to_str (fw)+") to \""+this->name+"\".\n"); f += fw; } } else { cvm::real const grad = this->dist2_lgrad (x, upper_wall); if (grad > 0.0) { fw = -0.5 * upper_wall_k * grad; if (cvm::debug()) cvm::log ("Applying an upper wall force ("+ cvm::to_str (fw)+") to \""+this->name+"\".\n"); f += fw; } } } if (tasks[task_Jacobian_force]) { cvm::increase_depth(); for (size_t i = 0; i < cvcs.size(); i++) { (cvcs[i])->calc_Jacobian_derivative(); } cvm::decrease_depth(); fj.reset(); for (size_t i = 0; i < cvcs.size(); i++) { // linear combination is assumed fj += 1.0 / ( cvm::real (cvcs.size()) * cvm::real ((cvcs[i])->sup_coeff) ) * (cvcs[i])->Jacobian_derivative(); } fj *= cvm::boltzmann() * cvm::temperature(); // the instantaneous Jacobian force was not included in the reported system force; // instead, it is subtracted from the applied force (silent Jacobian correction) if (! tasks[task_report_Jacobian_force]) f -= fj; } if (tasks[task_extended_lagrangian]) { cvm::real dt = cvm::dt(); // the total force is applied to the fictitious mass, while the // atoms only feel the harmonic force // fr: extended coordinate force; f: colvar force applied to atomic coordinates fr = f; fr += (-0.5 * ext_force_k) * this->dist2_lgrad (xr, x); f = (-0.5 * ext_force_k) * this->dist2_rgrad (xr, x); // leap frog: starting from x_i, f_i, v_(i-1/2) vr += (0.5 * dt) * fr / ext_mass; // Because of leapfrog, kinetic energy at time i is approximate kinetic_energy = 0.5 * ext_mass * vr * vr; potential_energy = 0.5 * ext_force_k * this->dist2(xr, x); // leap to v_(i+1/2) if (tasks[task_langevin]) { vr -= dt * ext_gamma * vr.real_value; vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass; } vr += (0.5 * dt) * fr / ext_mass; xr += dt * vr; xr.apply_constraints(); if (this->b_periodic) this->wrap (xr); } if (tasks[task_fdiff_velocity]) { // set it for the next step x_old = x; } if (cvm::debug()) cvm::log ("Done updating colvar \""+this->name+"\".\n"); return (potential_energy + kinetic_energy); } void colvar::communicate_forces() { if (cvm::debug()) cvm::log ("Communicating forces from colvar \""+this->name+"\".\n"); if (x.type() == colvarvalue::type_scalar) { for (size_t i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) * (std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1)) ); cvm::decrease_depth(); } } else { for (size_t i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff); cvm::decrease_depth(); } } if (cvm::debug()) cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n"); } // ******************** METRIC FUNCTIONS ******************** // Use the metrics defined by \link cvc \endlink objects bool colvar::periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) const { if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { cvm::fatal_error ("Error: requesting to histogram the " "collective variable \""+this->name+"\", but a " "pair of lower and upper boundaries must be " "defined.\n"); } if (period > 0.0) { if ( ((std::sqrt (this->dist2 (lb, ub))) / this->width) < 1.0E-10 ) { return true; } } return false; } bool colvar::periodic_boundaries() const { if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { cvm::fatal_error ("Error: requesting to histogram the " "collective variable \""+this->name+"\", but a " "pair of lower and upper boundaries must be " "defined.\n"); } return periodic_boundaries (lower_boundary, upper_boundary); } cvm::real colvar::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { return (cvcs[0])->dist2 (x1, x2); } colvarvalue colvar::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { return (cvcs[0])->dist2_lgrad (x1, x2); } colvarvalue colvar::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { return (cvcs[0])->dist2_rgrad (x1, x2); } cvm::real colvar::compare (colvarvalue const &x1, colvarvalue const &x2) const { return (cvcs[0])->compare (x1, x2); } void colvar::wrap (colvarvalue &x) const { (cvcs[0])->wrap (x); return; } // ******************** INPUT FUNCTIONS ******************** std::istream & colvar::read_restart (std::istream &is) { size_t const start_pos = is.tellg(); std::string conf; if ( !(is >> colvarparse::read_block ("colvar", conf)) ) { // this is not a colvar block is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } { std::string check_name = ""; if ( (get_keyval (conf, "name", check_name, std::string (""), colvarparse::parse_silent)) && (check_name != name) ) { cvm::fatal_error ("Error: the state file does not match the " "configuration file, at colvar \""+name+"\".\n"); } if (check_name.size() == 0) { cvm::fatal_error ("Error: Collective variable in the " "restart file without any identifier.\n"); } } if ( !(get_keyval (conf, "x", x, colvarvalue (x.type()), colvarparse::parse_silent)) ) { cvm::log ("Error: restart file does not contain " "the value of the colvar \""+ name+"\" .\n"); } else { cvm::log ("Restarting collective variable \""+name+"\" from value: "+ cvm::to_str (x)+"\n"); } if (tasks[colvar::task_extended_lagrangian]) { if ( !(get_keyval (conf, "extended_x", xr, colvarvalue (x.type()), colvarparse::parse_silent)) && !(get_keyval (conf, "extended_v", vr, colvarvalue (x.type()), colvarparse::parse_silent)) ) { cvm::log ("Error: restart file does not contain " "\"extended_x\" or \"extended_v\" for the colvar \""+ name+"\", but you requested \"extendedLagrangian\".\n"); } } if (tasks[task_extended_lagrangian]) { x_reported = xr; } else { x_reported = x; } if (tasks[task_output_velocity]) { if ( !(get_keyval (conf, "v", v_fdiff, colvarvalue (x.type()), colvarparse::parse_silent)) ) { cvm::log ("Error: restart file does not contain " "the velocity for the colvar \""+ name+"\", but you requested \"outputVelocity\".\n"); } if (tasks[task_extended_lagrangian]) { v_reported = vr; } else { v_reported = v_fdiff; } } return is; } std::istream & colvar::read_traj (std::istream &is) { size_t const start_pos = is.tellg(); if (tasks[task_output_value]) { if (!(is >> x)) { cvm::log ("Error: in reading the value of colvar \""+ this->name+"\" from trajectory.\n"); is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } if (tasks[task_extended_lagrangian]) { is >> xr; x_reported = xr; } else { x_reported = x; } } if (tasks[task_output_velocity]) { is >> v_fdiff; if (tasks[task_extended_lagrangian]) { is >> vr; v_reported = vr; } else { v_reported = v_fdiff; } } if (tasks[task_output_system_force]) { is >> ft; if (tasks[task_extended_lagrangian]) { is >> fr; ft_reported = fr; } else { ft_reported = ft; } } if (tasks[task_output_applied_force]) { is >> f; } return is; } // ******************** OUTPUT FUNCTIONS ******************** std::ostream & colvar::write_restart (std::ostream &os) { os << "colvar {\n" << " name " << name << "\n" << " x " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << x << "\n"; if (tasks[task_output_velocity]) { os << " v " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << v_reported << "\n"; } if (tasks[task_extended_lagrangian]) { os << " extended_x " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << xr << "\n" << " extended_v " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << vr << "\n"; } os << "}\n\n"; return os; } std::ostream & colvar::write_traj_label (std::ostream & os) { size_t const this_cv_width = x.output_width (cvm::cv_width); os << " "; if (tasks[task_output_value]) { os << " " << cvm::wrap_string (this->name, this_cv_width); if (tasks[task_extended_lagrangian]) { // restraint center os << " r_" << cvm::wrap_string (this->name, this_cv_width-2); } } if (tasks[task_output_velocity]) { os << " v_" << cvm::wrap_string (this->name, this_cv_width-2); if (tasks[task_extended_lagrangian]) { // restraint center os << " vr_" << cvm::wrap_string (this->name, this_cv_width-3); } } if (tasks[task_output_energy]) { os << " Ep_" << cvm::wrap_string (this->name, this_cv_width-3) << " Ek_" << cvm::wrap_string (this->name, this_cv_width-3); } if (tasks[task_output_system_force]) { os << " fs_" << cvm::wrap_string (this->name, this_cv_width-2); if (tasks[task_extended_lagrangian]) { // restraint center os << " fr_" << cvm::wrap_string (this->name, this_cv_width-3); } } if (tasks[task_output_applied_force]) { os << " fa_" << cvm::wrap_string (this->name, this_cv_width-2); } return os; } std::ostream & colvar::write_traj (std::ostream &os) { os << " "; if (tasks[task_output_value]) { if (tasks[task_extended_lagrangian]) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << x; } os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << x_reported; } if (tasks[task_output_velocity]) { if (tasks[task_extended_lagrangian]) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << v_fdiff; } os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << v_reported; } if (tasks[task_output_energy]) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << potential_energy << " " << kinetic_energy; } if (tasks[task_output_system_force]) { if (tasks[task_extended_lagrangian]) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << ft; } os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << ft_reported; } if (tasks[task_output_applied_force]) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << f; } return os; } void colvar::write_output_files() { if (cvm::b_analysis) { if (acf.size()) { cvm::log ("Writing acf to file \""+acf_outfile+"\".\n"); std::ofstream acf_os (acf_outfile.c_str()); if (! acf_os.good()) cvm::fatal_error ("Cannot open file \""+acf_outfile+"\".\n"); write_acf (acf_os); acf_os.close(); } if (runave_os.good()) { runave_os.close(); } } } // ******************** ANALYSIS FUNCTIONS ******************** void colvar::analyse() { if (tasks[task_runave]) { calc_runave(); } if (tasks[task_corrfunc]) { calc_acf(); } } inline void history_add_value (size_t const &history_length, std::list &history, colvarvalue const &new_value) { history.push_front (new_value); if (history.size() > history_length) history.pop_back(); } inline void history_incr (std::list< std::list > &history, std::list< std::list >::iterator &history_p) { if ((++history_p) == history.end()) history_p = history.begin(); } void colvar::calc_acf() { // using here an acf_stride-long list of vectors for either // coordinates (acf_x_history) or velocities (acf_v_history); each vector can // contain up to acf_length values, which are contiguous in memory // representation but separated by acf_stride in the time series; // the pointer to each vector is changed at every step if (! (acf_x_history.size() || acf_v_history.size()) ) { // first-step operations colvar *cfcv = (acf_colvar_name.size() ? cvm::colvar_p (acf_colvar_name) : this); if (cfcv->type() != this->type()) cvm::fatal_error ("Error: correlation function between \""+cfcv->name+ "\" and \""+this->name+"\" cannot be calculated, " "because their value types are different.\n"); acf_nframes = 0; cvm::log ("Colvar \""+this->name+"\": initializing ACF calculation.\n"); if (acf.size() < acf_length+1) acf.resize (acf_length+1, 0.0); switch (acf_type) { case acf_vel: // allocate space for the velocities history for (size_t i = 0; i < acf_stride; i++) { acf_v_history.push_back (std::list()); } acf_v_history_p = acf_v_history.begin(); break; case acf_coor: case acf_p2coor: // allocate space for the coordinates history for (size_t i = 0; i < acf_stride; i++) { acf_x_history.push_back (std::list()); } acf_x_history_p = acf_x_history.begin(); break; default: break; } } else { colvar *cfcv = (acf_colvar_name.size() ? cvm::colvar_p (acf_colvar_name) : this); switch (acf_type) { case acf_vel: if (tasks[task_fdiff_velocity]) { // calc() should do this already, but this only happens in a // simulation; better do it again in case a trajectory is // being read v_reported = v_fdiff = fdiff_velocity (x_old, cfcv->value()); } calc_vel_acf ((*acf_v_history_p), cfcv->velocity()); // store this value in the history history_add_value (acf_length+acf_offset, *acf_v_history_p, cfcv->velocity()); // if stride is larger than one, cycle among different histories history_incr (acf_v_history, acf_v_history_p); break; case acf_coor: calc_coor_acf ((*acf_x_history_p), cfcv->value()); history_add_value (acf_length+acf_offset, *acf_x_history_p, cfcv->value()); history_incr (acf_x_history, acf_x_history_p); break; case acf_p2coor: calc_p2coor_acf ((*acf_x_history_p), cfcv->value()); history_add_value (acf_length+acf_offset, *acf_x_history_p, cfcv->value()); history_incr (acf_x_history, acf_x_history_p); break; default: break; } } if (tasks[task_fdiff_velocity]) { // set it for the next step x_old = x; } } void colvar::calc_vel_acf (std::list &v_list, colvarvalue const &v) { // loop over stored velocities and add to the ACF, but only the // length is sufficient to hold an entire row of ACF values if (v_list.size() >= acf_length+acf_offset) { std::list::iterator vs_i = v_list.begin(); std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) vs_i++; // current vel with itself *(acf_i++) += v.norm2(); // inner products of previous velocities with current (acf_i and // vs_i are updated) colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i); acf_nframes++; } } void colvar::calc_coor_acf (std::list &x_list, colvarvalue const &x) { // same as above but for coordinates if (x_list.size() >= acf_length+acf_offset) { std::list::iterator xs_i = x_list.begin(); std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) xs_i++; *(acf_i++) += x.norm2(); colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i); acf_nframes++; } } void colvar::calc_p2coor_acf (std::list &x_list, colvarvalue const &x) { // same as above but with second order Legendre polynomial instead // of just the scalar product if (x_list.size() >= acf_length+acf_offset) { std::list::iterator xs_i = x_list.begin(); std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) xs_i++; // value of P2(0) = 1 *(acf_i++) += 1.0; colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i); acf_nframes++; } } void colvar::write_acf (std::ostream &os) { if (!acf_nframes) cvm::log ("Warning: ACF was not calculated (insufficient frames).\n"); os.setf (std::ios::scientific, std::ios::floatfield); os << "# Autocorrelation function for collective variable \"" << this->name << "\"\n"; // one frame is used for normalization, the statistical sample is // hence decreased os << "# nframes = " << (acf_normalize ? acf_nframes - 1 : acf_nframes) << "\n"; cvm::real const acf_norm = acf.front() / cvm::real (acf_nframes); std::vector::iterator acf_i; size_t it = acf_offset; for (acf_i = acf.begin(); acf_i != acf.end(); acf_i++) { os << std::setw (cvm::it_width) << acf_stride * (it++) << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << ( acf_normalize ? (*acf_i)/(acf_norm * cvm::real (acf_nframes)) : (*acf_i)/(cvm::real (acf_nframes)) ) << "\n"; } } void colvar::calc_runave() { if (!x_history.size()) { runave.type (x.type()); runave.reset(); // first-step operations if (cvm::debug()) cvm::log ("Colvar \""+this->name+ "\": initializing running average calculation.\n"); acf_nframes = 0; x_history.push_back (std::list()); x_history_p = x_history.begin(); } else { if ( (cvm::step_relative() % runave_stride) == 0) { if ((*x_history_p).size() >= runave_length-1) { runave = x; for (std::list::iterator xs_i = (*x_history_p).begin(); xs_i != (*x_history_p).end(); xs_i++) { runave += (*xs_i); } runave *= 1.0 / cvm::real (runave_length); runave.apply_constraints(); runave_variance = 0.0; runave_variance += this->dist2 (x, runave); for (std::list::iterator xs_i = (*x_history_p).begin(); xs_i != (*x_history_p).end(); xs_i++) { runave_variance += this->dist2 (x, (*xs_i)); } runave_variance *= 1.0 / cvm::real (runave_length-1); runave_os << std::setw (cvm::it_width) << cvm::step_relative() << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << runave << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << std::sqrt (runave_variance) << "\n"; } history_add_value (runave_length, *x_history_p, x); } } } diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 7a3ffa190..1588fcee8 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -1,781 +1,807 @@ #include "colvarmodule.h" #include "colvarparse.h" #include "colvaratoms.h" -// atom member functions depend tightly on the MD interface, and are +// member functions of the "atom" class depend tightly on the MD interface, and are // thus defined in colvarproxy_xxx.cpp +// in this file only atom_group functions are defined -// atom_group member functions +// Note: "conf" is the configuration of the cvc who is using this atom group; +// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...) cvm::atom_group::atom_group (std::string const &conf, - char const *key, - atom_group *ref_pos_group_in) + char const *key) : b_center (false), b_rotate (false), b_user_defined_fit (false), - ref_pos_group (NULL), // this is always set within parse(), - // regardless of ref_pos_group_in + ref_pos_group (NULL), + b_fit_gradients (false), noforce (false) { cvm::log ("Defining atom group \""+ std::string (key)+"\".\n"); - parse (conf, key, ref_pos_group_in); + // real work is done by parse + parse (conf, key); +} + + +cvm::atom_group::atom_group (std::vector const &atoms) + : b_dummy (false), b_center (false), b_rotate (false), + ref_pos_group (NULL), b_fit_gradients (false), + noforce (false) +{ + this->reserve (atoms.size()); + for (size_t i = 0; i < atoms.size(); i++) { + this->push_back (atoms[i]); + } + total_mass = 0.0; + for (cvm::atom_iter ai = this->begin(); + ai != this->end(); ai++) { + total_mass += ai->mass; + } +} + + +cvm::atom_group::atom_group() + : b_dummy (false), b_center (false), b_rotate (false), + ref_pos_group (NULL), b_fit_gradients (false), + noforce (false) +{ + total_mass = 0.0; +} + + +cvm::atom_group::~atom_group() +{ + if (ref_pos_group) { + delete ref_pos_group; + ref_pos_group = NULL; + } +} + + +void cvm::atom_group::add_atom (cvm::atom const &a) +{ + if (b_dummy) { + cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n"); + } else { + this->push_back (a); + total_mass += a.mass; + } } void cvm::atom_group::parse (std::string const &conf, - char const *key, - atom_group *ref_pos_group_in) + char const *key) { std::string group_conf; // save_delimiters is set to false for this call, because "conf" is // not the config string of this group, but of its parent object // (which has already taken care of the delimiters) save_delimiters = false; key_lookup (conf, key, group_conf, dummy_pos); // restoring the normal value, because we do want keywords checked // inside "group_conf" save_delimiters = true; if (group_conf.size() == 0) { cvm::fatal_error ("Error: atom group \""+ std::string (key)+"\" is set, but " "has no definition.\n"); } cvm::increase_depth(); cvm::log ("Initializing atom group \""+std::string (key)+"\".\n"); // whether or not to include messages in the log // colvarparse::Parse_Mode mode = parse_silent; // { // bool b_verbose; // get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent); // if (b_verbose) mode = parse_normal; // } colvarparse::Parse_Mode mode = parse_normal; { // get the atoms by numbers std::vector atom_indexes; if (get_keyval (group_conf, "atomNumbers", atom_indexes, atom_indexes, mode)) { if (atom_indexes.size()) { this->reserve (this->size()+atom_indexes.size()); for (size_t i = 0; i < atom_indexes.size(); i++) { this->push_back (cvm::atom (atom_indexes[i])); } } else cvm::fatal_error ("Error: no numbers provided for \"" "atomNumbers\".\n"); } } { std::string range_conf = ""; size_t pos = 0; while (key_lookup (group_conf, "atomNumbersRange", range_conf, pos)) { if (range_conf.size()) { std::istringstream is (range_conf); int initial, final; char dash; if ( (is >> initial) && (initial > 0) && (is >> dash) && (dash == '-') && (is >> final) && (final > 0) ) { for (int anum = initial; anum <= final; anum++) { this->push_back (cvm::atom (anum)); } range_conf = ""; continue; } } cvm::fatal_error ("Error: no valid definition for \"" "atomNumbersRange\", \""+ range_conf+"\".\n"); } } std::vector psf_segids; get_keyval (group_conf, "psfSegID", psf_segids, std::vector (), mode); for (std::vector::iterator psii = psf_segids.begin(); psii < psf_segids.end(); psii++) { if ( (psii->size() == 0) || (psii->size() > 4) ) { cvm::fatal_error ("Error: invalid segmend identifier provided, \""+ (*psii)+"\".\n"); } } { std::string range_conf = ""; size_t pos = 0; size_t range_count = 0; std::vector::iterator psii = psf_segids.begin(); while (key_lookup (group_conf, "atomNameResidueRange", range_conf, pos)) { range_count++; if (range_count > psf_segids.size()) { cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than " "values of \"psfSegID\".\n"); } std::string const &psf_segid = psf_segids.size() ? *psii : std::string (""); if (range_conf.size()) { std::istringstream is (range_conf); std::string atom_name; int initial, final; char dash; if ( (is >> atom_name) && (atom_name.size()) && (is >> initial) && (initial > 0) && (is >> dash) && (dash == '-') && (is >> final) && (final > 0) ) { for (int resid = initial; resid <= final; resid++) { this->push_back (cvm::atom (resid, atom_name, psf_segid)); } range_conf = ""; } else { cvm::fatal_error ("Error: cannot parse definition for \"" "atomNameResidueRange\", \""+ range_conf+"\".\n"); } } else { cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n"); } if (psf_segid.size()) psii++; } } { // read the atoms from a file std::string atoms_file_name; if (get_keyval (group_conf, "atomsFile", atoms_file_name, std::string (""), mode)) { std::string atoms_col; if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) { cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n"); } double atoms_col_value; bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode); if (atoms_col_value_defined && (!atoms_col_value)) cvm::fatal_error ("Error: atomsColValue, " "if provided, must be non-zero.\n"); cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value); } } for (std::vector::iterator a1 = this->begin(); a1 != this->end(); a1++) { std::vector::iterator a2 = a1; ++a2; for ( ; a2 != this->end(); a2++) { if (a1->id == a2->id) { if (cvm::debug()) cvm::log ("Discarding doubly counted atom with number "+ cvm::to_str (a1->id+1)+".\n"); a2 = this->erase (a2); if (a2 == this->end()) break; } } } if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) { b_dummy = true; this->total_mass = 1.0; } else b_dummy = false; if (b_dummy && (this->size())) cvm::fatal_error ("Error: cannot set up group \""+ std::string (key)+"\" as a dummy atom " "and provide it with atom definitions.\n"); #if (! defined (COLVARS_STANDALONE)) if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) { cvm::fatal_error ("Error: no atoms defined for atom group \""+ std::string (key)+"\".\n"); } #endif if (!b_dummy) { this->total_mass = 0.0; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { this->total_mass += ai->mass; } } if (!b_dummy) get_keyval (group_conf, "disableForces", noforce, false, mode); // FITTING OPTIONS bool b_defined_center = get_keyval (group_conf, "centerReference", b_center, false, mode); bool b_defined_rotate = get_keyval (group_conf, "rotateReference", b_rotate, false, mode); + b_fit_gradients = get_keyval (group_conf, "fitGradients", b_fit_gradients, false, mode); + // this cannot be shortened to one statement because lazy evaluation may prevent one // function from being called! b_user_defined_fit = b_defined_center || b_defined_rotate; // if ((b_center || b_rotate) && b_dummy) // cvm::fatal_error ("Error: cannot set \"centerReference\" or " // "\"rotateReference\" when \"dummyAtom\" is defined.\n"); if (b_center || b_rotate) { - // instead of this group, define another group (refPositionsGroup) to be the one - // used to fit the coordinates if (key_lookup (group_conf, "refPositionsGroup")) { + // instead of this group, define another group (refPositionsGroup) + // to be the one used to fit the coordinates if (ref_pos_group) { cvm::fatal_error ("Error: the atom group \""+ std::string (key)+"\" has already a reference group " "for the rototranslational fit, which was communicated by the " "colvar component. You should not use refPositionsGroup " "in this case.\n"); } cvm::log ("Within atom group \""+std::string (key)+"\":\n"); ref_pos_group = new atom_group (group_conf, "refPositionsGroup"); } atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; - if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) { - if (ref_pos.size() != group_for_fit->size()) { - cvm::fatal_error ("Error: the number of reference positions provided ("+ - cvm::to_str (ref_pos.size())+ - ") does not match the number of atoms of group \""+ - std::string (key)+ - "\" ("+cvm::to_str (group_for_fit->size())+").\n"); - } - } + get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode); std::string ref_pos_file; if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) { if (ref_pos.size()) { cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and " "\"refPositions\" at the same time.\n"); } std::string ref_pos_col; double ref_pos_col_value; if (get_keyval (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) { // if provided, use PDB column to select coordinates bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode); if (found && !ref_pos_col_value) cvm::fatal_error ("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); } else { // if not, rely on existing atom indices for the group group_for_fit->create_sorted_ids(); } cvm::load_coords (ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids, ref_pos_col, ref_pos_col_value); } if (ref_pos.size()) { if (b_rotate) { if (ref_pos.size() != group_for_fit->size()) cvm::fatal_error ("Error: the number of reference positions provided ("+ cvm::to_str (ref_pos.size())+ ") does not match the number of atoms within \""+ std::string (key)+ "\" ("+cvm::to_str (group_for_fit->size())+ "): to perform a rotational fit, "+ "these numbers should be equal.\n"); } // save the center of geometry of ref_pos and subtract it center_ref_pos(); } else { #if (! defined (COLVARS_STANDALONE)) cvm::fatal_error ("Error: no reference positions provided.\n"); #endif } + if (b_fit_gradients) { + group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::atom_pos (0.0, 0.0, 0.0)); + rot.request_group1_gradients (group_for_fit->size()); + } + if (b_rotate && !noforce) { cvm::log ("Warning: atom group \""+std::string (key)+ - "\" will be fitted automatically onto a fixed orientation: " - "in few cases, torques applied on this group " - "may make the simulation unstable. " - "If this happens, set \"disableForces\" to yes " - "for this group.\n"); + "\" will be aligned to a fixed orientation given by the reference positions provided. " + "If the internal structure of the group changes too much (i.e. its RMSD is comparable " + "to its radius of gyration), the optimal rotation and its gradients may become discontinuous. " + "If that happens, use refPositionsGroup (or a different definition for it if already defined) " + "to align the coordinates.\n"); // initialize rot member data rot.request_group1_gradients (this->size()); } + } if (cvm::debug()) cvm::log ("Done initializing atom group with name \""+ std::string (key)+"\".\n"); this->check_keywords (group_conf, key); cvm::log ("Atom group \""+std::string (key)+"\" defined, "+ cvm::to_str (this->size())+" initialized: total mass = "+ cvm::to_str (this->total_mass)+".\n"); cvm::decrease_depth(); } -cvm::atom_group::atom_group (std::vector const &atoms) - : b_dummy (false), b_center (false), b_rotate (false), - ref_pos_group (NULL), noforce (false) -{ - this->reserve (atoms.size()); - for (size_t i = 0; i < atoms.size(); i++) { - this->push_back (atoms[i]); - } - total_mass = 0.0; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - total_mass += ai->mass; - } -} - - -cvm::atom_group::atom_group() - : b_dummy (false), b_center (false), b_rotate (false), - ref_pos_group (NULL), noforce (false) -{ - total_mass = 0.0; -} - - -cvm::atom_group::~atom_group() -{ - if (ref_pos_group) { - delete ref_pos_group; - ref_pos_group = NULL; - } -} - - -void cvm::atom_group::add_atom (cvm::atom const &a) -{ - if (b_dummy) { - cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n"); - } else { - this->push_back (a); - total_mass += a.mass; - } -} - - void cvm::atom_group::create_sorted_ids (void) { // Only do the work if the vector is not yet populated if (sorted_ids.size()) return; std::list temp_id_list; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { temp_id_list.push_back (ai->id); } temp_id_list.sort(); temp_id_list.unique(); if (temp_id_list.size() != this->size()) { cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " + cvm::to_str(temp_id_list.size()) + " unique atom IDs instead of" + cvm::to_str(this->size()) + ").\n"); } sorted_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); return; } void cvm::atom_group::center_ref_pos() { // save the center of geometry of ref_pos and then subtract it from // them; in this way it will be possible to use ref_pos also for // the rotational fit // This is called either by atom_group::parse or by CVCs that set // reference positions (eg. RMSD, eigenvector) ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0); std::vector::iterator pi = ref_pos.begin(); for ( ; pi != ref_pos.end(); pi++) { ref_pos_cog += *pi; } ref_pos_cog /= (cvm::real) ref_pos.size(); for (std::vector::iterator pi = ref_pos.begin(); pi != ref_pos.end(); pi++) { (*pi) -= ref_pos_cog; } } void cvm::atom_group::read_positions() { if (b_dummy) return; -#if (! defined (COLVARS_STANDALONE)) - if (!this->size()) - cvm::fatal_error ("Error: no atoms defined in the requested group.\n"); -#endif - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_position(); } if (ref_pos_group) ref_pos_group->read_positions(); +} +void cvm::atom_group::calc_apply_roto_translation() +{ atom_group *fit_group = ref_pos_group ? ref_pos_group : this; if (b_center) { - // store aside the current center of geometry (all positions will be - // set to the closest images to the first one) and then center on - // the origin + // center on the origin first cvm::atom_pos const cog = fit_group->center_of_geometry(); for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos -= cog; } } if (b_rotate) { // rotate the group (around the center of geometry if b_center is - // true, around the origin otherwise); store the rotation, in - // order to bring back the forces to the original frame before - // applying them + // true, around the origin otherwise) rot.calc_optimal_rotation (fit_group->positions(), ref_pos); for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos = rot.rotate (ai->pos); } } if (b_center) { - // use the center of geometry of ref_pos + // align with the center of geometry of ref_pos for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos += ref_pos_cog; } } } - void cvm::atom_group::apply_translation (cvm::rvector const &t) { if (b_dummy) return; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos += t; } } - void cvm::atom_group::apply_rotation (cvm::rotation const &rot) { if (b_dummy) return; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos = rot.rotate (ai->pos); } } void cvm::atom_group::read_velocities() { if (b_dummy) return; if (b_rotate) { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_velocity(); ai->vel = rot.rotate (ai->vel); } } else { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_velocity(); } } } void cvm::atom_group::read_system_forces() { if (b_dummy) return; if (b_rotate) { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_system_force(); ai->system_force = rot.rotate (ai->system_force); } } else { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_system_force(); } } } -/* This is deprecated. -cvm::atom_pos cvm::atom_group::center_of_geometry (cvm::atom_pos const &ref_pos) -{ - if (b_dummy) { - cvm::select_closest_image (dummy_atom_pos, ref_pos); - return dummy_atom_pos; - } - - cvm::atom_pos cog (0.0, 0.0, 0.0); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - cvm::select_closest_image (ai->pos, ref_pos); - cog += ai->pos; - } - cog /= this->size(); - return cog; -} */ - cvm::atom_pos cvm::atom_group::center_of_geometry() const { if (b_dummy) return dummy_atom_pos; cvm::atom_pos cog (0.0, 0.0, 0.0); for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { cog += ai->pos; } cog /= this->size(); return cog; } -/* This is deprecated. -cvm::atom_pos cvm::atom_group::center_of_mass (cvm::atom_pos const &ref_pos) -{ - if (b_dummy) { - cvm::select_closest_image (dummy_atom_pos, ref_pos); - return dummy_atom_pos; - } - - cvm::atom_pos com (0.0, 0.0, 0.0); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - cvm::select_closest_image (ai->pos, ref_pos); - com += ai->mass * ai->pos; - } - com /= this->total_mass; - return com; -}*/ - cvm::atom_pos cvm::atom_group::center_of_mass() const { if (b_dummy) return dummy_atom_pos; cvm::atom_pos com (0.0, 0.0, 0.0); for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { com += ai->mass * ai->pos; } com /= this->total_mass; return com; } void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad) { if (b_dummy) return; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->grad = (ai->mass/this->total_mass) * grad; } } +void cvm::atom_group::calc_fit_gradients() +{ + if (b_dummy) return; + + if (cvm::debug()) + cvm::log ("Calculating fit gradients.\n"); + + atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; + group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::rvector (0.0, 0.0, 0.0)); + + if (b_center) { + // add the center of geometry contribution to the gradients + for (size_t i = 0; i < this->size(); i++) { + // need to bring the gradients in original frame first + cvm::rvector const atom_grad = b_rotate ? + (rot.inverse()).rotate ((*this)[i].grad) : + (*this)[i].grad; + for (size_t j = 0; j < group_for_fit->size(); j++) { + group_for_fit->fit_gradients[j] += + (-1.0)/(cvm::real (group_for_fit->size())) * + atom_grad; + } + } + } + + if (b_rotate) { + + // add the rotation matrix contribution to the gradients + cvm::rotation const rot_inv = rot.inverse(); + cvm::atom_pos const cog = this->center_of_geometry(); + + for (size_t i = 0; i < this->size(); i++) { + + cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos))); + + for (size_t j = 0; j < group_for_fit->size(); j++) { + // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i + cvm::quaternion const dxdq = + rot.q.position_derivative_inner (pos_orig, (*this)[i].grad); + // multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients + for (size_t iq = 0; iq < 4; iq++) { + group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq]; + } + } + } + } + if (cvm::debug()) + cvm::log ("Done calculating fit gradients.\n"); +} + + std::vector cvm::atom_group::positions() const { if (b_dummy) cvm::fatal_error ("Error: positions are not available " "from a dummy atom group.\n"); std::vector x (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator xi = x.begin(); for ( ; ai != this->end(); xi++, ai++) { *xi = ai->pos; } return x; } std::vector cvm::atom_group::positions_shifted (cvm::rvector const &shift) const { if (b_dummy) cvm::fatal_error ("Error: positions are not available " "from a dummy atom group.\n"); std::vector x (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator xi = x.begin(); for ( ; ai != this->end(); xi++, ai++) { *xi = (ai->pos + shift); } return x; } std::vector cvm::atom_group::velocities() const { if (b_dummy) cvm::fatal_error ("Error: velocities are not available " "from a dummy atom group.\n"); std::vector v (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator vi = v.begin(); for ( ; ai != this->end(); vi++, ai++) { *vi = ai->vel; } return v; } std::vector cvm::atom_group::system_forces() const { if (b_dummy) cvm::fatal_error ("Error: system forces are not available " "from a dummy atom group.\n"); std::vector f (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator fi = f.begin(); for ( ; ai != this->end(); fi++, ai++) { *fi = ai->system_force; } return f; } cvm::rvector cvm::atom_group::system_force() const { if (b_dummy) cvm::fatal_error ("Error: system forces are not available " "from a dummy atom group.\n"); cvm::rvector f (0.0); for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { f += ai->system_force; } return f; } void cvm::atom_group::apply_colvar_force (cvm::real const &force) { if (b_dummy) return; if (noforce) cvm::fatal_error ("Error: sending a force to a group that has " "\"disableForces\" defined.\n"); if (b_rotate) { - // get the forces back from the rotated frame + // rotate forces back to the original frame cvm::rotation const rot_inv = rot.inverse(); for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force (rot_inv.rotate (force * ai->grad)); } } else { - // no need to manipulate gradients, they are still in the original - // frame for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force (force * ai->grad); } } + + if (b_fit_gradients) { + + atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; + + // add the contribution from the roto-translational fit to the gradients + if (b_rotate) { + // rotate forces back to the original frame + cvm::rotation const rot_inv = rot.inverse(); + for (size_t j = 0; j < group_for_fit->size(); j++) { + (*group_for_fit)[j].apply_force (rot_inv.rotate (force * fit_gradients[j])); + } + } else { + for (size_t j = 0; j < group_for_fit->size(); j++) { + (*group_for_fit)[j].apply_force (force * fit_gradients[j]); + } + } + } + } void cvm::atom_group::apply_force (cvm::rvector const &force) { if (b_dummy) return; if (noforce) cvm::fatal_error ("Error: sending a force to a group that has " "\"disableForces\" defined.\n"); if (b_rotate) { cvm::rotation const rot_inv = rot.inverse(); for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force (rot_inv.rotate ((ai->mass/this->total_mass) * force)); } } else { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force ((ai->mass/this->total_mass) * force); } } } void cvm::atom_group::apply_forces (std::vector const &forces) { if (b_dummy) return; if (noforce) cvm::fatal_error ("Error: sending a force to a group that has " "\"disableForces\" defined.\n"); if (forces.size() != this->size()) { cvm::fatal_error ("Error: trying to apply an array of forces to an atom " "group which does not have the same length.\n"); } if (b_rotate) { cvm::rotation const rot_inv = rot.inverse(); cvm::atom_iter ai = this->begin(); std::vector::const_iterator fi = forces.begin(); for ( ; ai != this->end(); fi++, ai++) { ai->apply_force (rot_inv.rotate (*fi)); } } else { cvm::atom_iter ai = this->begin(); std::vector::const_iterator fi = forces.begin(); for ( ; ai != this->end(); fi++, ai++) { ai->apply_force (*fi); } } } - diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 6ffda6db5..d2be1dfe6 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -1,322 +1,327 @@ // -*- c++ -*- #ifndef COLVARATOMS_H #define COLVARATOMS_H #include "colvarmodule.h" #include "colvarparse.h" /// \brief Stores numeric id, mass and all mutable data for an atom, /// mostly used by a \link cvc \endlink /// /// This class may be used (although not necessarily) to keep atomic /// data (id, mass, position and collective variable derivatives) /// altogether. There may be multiple instances with identical /// numeric id, all acting independently: forces communicated through /// these instances will be summed together. /// /// Read/write operations depend on the underlying code: hence, some /// member functions are defined in colvarproxy_xxx.h. class colvarmodule::atom { protected: /// \brief Index in the list of atoms involved by the colvars (\b /// NOT in the global topology!) size_t index; public: /// Internal identifier (zero-based) int id; /// Mass cvm::real mass; /// \brief Current position (copied from the program, can be /// manipulated) cvm::atom_pos pos; /// \brief Current velocity (copied from the program, can be /// manipulated) cvm::rvector vel; /// \brief System force at the previous step (copied from the /// program, can be manipulated) cvm::rvector system_force; /// \brief Gradient of a scalar collective variable with respect /// to this atom /// /// This can only handle a scalar collective variable (i.e. when /// the \link colvarvalue::real_value \endlink member is used /// from the \link colvarvalue \endlink class), which is also the /// most frequent case. For more complex types of \link /// colvarvalue \endlink objects, atomic gradients should be /// defined within the specific \link cvc \endlink /// implementation cvm::rvector grad; /// \brief Default constructor, setting id to a non-valid value inline atom() {} /// \brief Initialize an atom for collective variable calculation /// and get its internal identifier \param atom_number Atom index in /// the system topology (starting from 1) atom (int const &atom_number); /// \brief Initialize an atom for collective variable calculation /// and get its internal identifier \param residue Residue number /// \param atom_name Name of the atom in the residue \param /// segment_id For PSF topologies, the segment identifier; for other /// type of topologies, may not be required atom (cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id = std::string ("")); /// Copy constructor atom (atom const &a); /// Destructor ~atom(); /// Set non-constant data (everything except id and mass) to zero inline void reset_data() { pos = atom_pos (0.0); vel = grad = system_force = rvector (0.0); } /// Get the current position void read_position(); /// Get the current velocity void read_velocity(); /// Get the system force void read_system_force(); /// \brief Apply a force to the atom /// /// The force will be used later by the MD integrator, the /// collective variables module does not integrate equations of /// motion. Multiple calls to this function by either the same /// \link atom \endlink object or different objects with identical /// \link id \endlink, will all add to the existing MD force. void apply_force (cvm::rvector const &new_force); }; /// \brief Group of \link atom \endlink objects, mostly used by a /// \link cvc \endlink /// /// This class inherits from \link colvarparse \endlink and from /// std::vector, and hence all functions and /// operators (including the bracket operator, group[i]) can be used /// on an \link atom_group \endlink object. It can be initialized as /// a vector, or by parsing a keyword in the configuration. class colvarmodule::atom_group : public std::vector, public colvarparse { public: // Note: all members here are kept public, to make possible to any // object accessing and manipulating them /// \brief If this option is on, this group merely acts as a wrapper /// for a fixed position; any calls to atoms within or to /// functions that return disaggregated data will fail bool b_dummy; /// \brief dummy atom position cvm::atom_pos dummy_atom_pos; /// Sorted list of zero-based (internal) atom ids /// (populated on-demand by create_sorted_ids) std::vector sorted_ids; /// Allocates and populates the sorted list of atom ids void create_sorted_ids (void); - /// \brief Before calculating colvars, move the group to overlap the - /// center of mass of reference coordinates + /// \brief When updating atomic coordinates, translate them to align with the + /// center of mass of the reference coordinates bool b_center; - /// \brief Right after updating atom coordinates (and after - /// centering coordinates, if b_center is true), rotate the group to - /// overlap the reference coordinates. You should not manipulate - /// atoms individually if you turn on this flag. + /// \brief When updating atom coordinates (and after + /// centering them if b_center is set), rotate the group to + /// align with the reference coordinates. /// /// Note: gradients will be calculated in the rotated frame: when /// forces will be applied, they will rotated back to the original /// frame bool b_rotate; - /// Rotation between the group and its reference coordinates + /// The rotation calculated automatically if b_rotate is defined cvm::rotation rot; - /// \brief Indicates that the user has specified centerReference or - /// rotateReference (and the related reference coordinates): - /// cvc's (eg rmsd, eigenvector) should not override these settings + /// \brief Indicates that the user has explicitly set centerReference or + /// rotateReference, and the corresponding reference: + /// cvc's (eg rmsd, eigenvector) will not override the user's choice bool b_user_defined_fit; - /// \brief use these reference coordinates, for b_center, b_rotate, or to compute - /// certain colvar components (orientation, rmsd, etc) + /// \brief Whether or not the derivatives of the roto-translation + /// should be included when calculating the colvar's gradients (default: no) + bool b_fit_gradients; + + /// \brief use reference coordinates for b_center or b_rotate std::vector ref_pos; + /// \brief Center of geometry of the reference coordinates; regardless /// of whether b_center is true, ref_pos is centered to zero at /// initialization, and ref_pos_cog serves to center the positions cvm::atom_pos ref_pos_cog; - /// \brief In case b_center or b_rotate is true, fit this group to - /// the reference positions (default: the parent group itself) - atom_group *ref_pos_group; - + /// \brief If b_center or b_rotate is true, use this group to + /// define the transformation (default: this group itself) + atom_group *ref_pos_group; + /// Total mass of the atom group cvm::real total_mass; /// \brief Don't apply any force on this group (use its coordinates /// only to calculate a colvar) bool noforce; /// \brief Initialize the group by looking up its configuration /// string in conf and parsing it; this is actually done by parse(), /// which is a member function so that a group can be initialized /// also after construction atom_group (std::string const &conf, - char const *key, - atom_group *ref_pos_group = NULL); + char const *key); /// \brief Initialize the group by looking up its configuration /// string in conf and parsing it void parse (std::string const &conf, - char const *key, - atom_group *ref_pos_group = NULL); + char const *key); /// \brief Initialize the group after a temporary vector of atoms atom_group (std::vector const &atoms); /// \brief Add an atom to this group void add_atom (cvm::atom const &a); /// \brief Default constructor atom_group(); /// \brief Destructor ~atom_group(); /// \brief Get the current positions; if b_center or b_rotate are - /// true, center and/or rotate the coordinates right after reading - /// them + /// true, calc_apply_roto_translation() will be called too void read_positions(); + /// \brief (Re)calculate the optimal roto-translation + void calc_apply_roto_translation(); + /// \brief Save center of geometry fo ref positions, /// then subtract it void center_ref_pos(); /// \brief Move all positions void apply_translation (cvm::rvector const &t); /// \brief Rotate all positions void apply_rotation (cvm::rotation const &q); /// \brief Get the current velocities; this must be called always /// *after* read_positions(); if b_rotate is defined, the same /// rotation applied to the coordinates will be used void read_velocities(); /// \brief Get the current system_forces; this must be called always /// *after* read_positions(); if b_rotate is defined, the same /// rotation applied to the coordinates will be used void read_system_forces(); /// Call reset_data() for each atom inline void reset_atoms_data() { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) ai->reset_data(); + if (ref_pos_group) + ref_pos_group->reset_atoms_data(); } /// \brief Return a copy of the current atom positions std::vector positions() const; /// \brief Return a copy of the current atom positions, shifted by a constant vector std::vector positions_shifted (cvm::rvector const &shift) const; - /// \brief Return the center of geometry of the positions \param ref_pos - /// Use the closest periodic images to this position - cvm::atom_pos center_of_geometry (cvm::atom_pos const &ref_pos); /// \brief Return the center of geometry of the positions, assuming /// that coordinates are already pbc-wrapped cvm::atom_pos center_of_geometry() const; - /// \brief Return the center of mass of the positions \param ref_pos - /// Use the closest periodic images to this position - cvm::atom_pos center_of_mass (cvm::atom_pos const &ref_pos); /// \brief Return the center of mass of the positions, assuming that /// coordinates are already pbc-wrapped cvm::atom_pos center_of_mass() const; - /// \brief Store atom positions from the previous step + /// \brief Atom positions at the previous step std::vector old_pos; /// \brief Return a copy of the current atom velocities std::vector velocities() const; /// \brief Return a copy of the system forces std::vector system_forces() const; /// \brief Return a copy of the aggregated total force on the group cvm::rvector system_force() const; /// \brief Shorthand: save the specified gradient on each atom, /// weighting with the atom mass (mostly used in combination with /// \link center_of_mass() \endlink) void set_weighted_gradient (cvm::rvector const &grad); + /// \brief Calculate the derivatives of the fitting transformation + void calc_fit_gradients(); + + /// \brief Derivatives of the fitting transformation + std::vector fit_gradients; + /// \brief Used by a (scalar) colvar to apply its force on its \link /// atom_group \endlink members /// /// The (scalar) force is multiplied by the colvar gradient for each /// atom; this should be used when a colvar with scalar \link /// colvarvalue \endlink type is used (this is the most frequent /// case: for colvars with a non-scalar type, the most convenient /// solution is to sum together the Cartesian forces from all the /// colvar components, and use apply_force() or apply_forces()). If /// the group is being rotated to a reference frame (e.g. to express /// the colvar independently from the solute rotation), the - /// gradients are temporarily to the original frame. + /// gradients are temporarily rotated to the original frame. void apply_colvar_force (cvm::real const &force); /// \brief Apply a force "to the center of mass", i.e. the force is /// distributed on each atom according to its mass /// /// If the group is being rotated to a reference frame (e.g. to /// express the colvar independently from the solute rotation), the /// force is rotated back to the original frame. Colvar gradients /// are not used, either because they were not defined (e.g because /// the colvar has not a scalar value) or the biases require to /// micromanage the force. void apply_force (cvm::rvector const &force); /// \brief Apply an array of forces directly on the individual /// atoms; the length of the specified vector must be the same of /// this \link atom_group \endlink. /// /// If the group is being rotated to a reference frame (e.g. to /// express the colvar independently from the solute rotation), the /// forces are rotated back to the original frame. Colvar gradients /// are not used, either because they were not defined (e.g because /// the colvar has not a scalar value) or the biases require to /// micromanage the forces. void apply_forces (std::vector const &forces); }; #endif diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index a3a483d60..2d03b5c8d 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -1,1709 +1,1713 @@ #include #include #include #include #include #include // used to set the absolute path of a replica file #if defined (WIN32) && !defined(__CYGWIN__) #include #define CHDIR ::_chdir #define GETCWD ::_getcwd #define PATHSEP "\\" #else #include #define CHDIR ::chdir #define GETCWD ::getcwd #define PATHSEP "/" #endif #include "colvar.h" #include "colvarbias_meta.h" colvarbias_meta::colvarbias_meta() : colvarbias(), state_file_step (0), new_hills_begin (hills.end()) { } colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) : colvarbias (conf, key), state_file_step (0), new_hills_begin (hills.end()) { if (cvm::n_abf_biases > 0) cvm::log ("Warning: ABF and metadynamics running at the " "same time can lead to inconsistent results.\n"); get_keyval (conf, "hillWeight", hill_weight, 0.01); if (hill_weight == 0.0) cvm::log ("Warning: hillWeight has been set to zero, " "this bias will have no effect.\n"); get_keyval (conf, "newHillFrequency", new_hill_freq, 1000); get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0); + { + bool b_replicas = false; + get_keyval (conf, "multipleReplicas", b_replicas, false); + if (b_replicas) + comm = multiple_replicas; + else + comm = single_replica; + } + get_keyval (conf, "useGrids", use_grids, true); if (use_grids) { get_keyval (conf, "gridsUpdateFrequency", grids_freq, new_hill_freq); get_keyval (conf, "rebinGrids", rebin_grids, false); expand_grids = false; for (size_t i = 0; i < colvars.size(); i++) { if (colvars[i]->expand_boundaries) { expand_grids = true; cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": Will expand grids when the colvar \""+ colvars[i]->name+"\" approaches its boundaries.\n"); } } get_keyval (conf, "keepHills", keep_hills, false); get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true); get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false); for (size_t i = 0; i < colvars.size(); i++) { colvars[i]->enable (colvar::task_grid); } hills_energy = new colvar_grid_scalar (colvars); hills_energy_gradients = new colvar_grid_gradient (colvars); } else { rebin_grids = false; keep_hills = false; dump_fes = false; dump_fes_save = false; dump_replica_fes = false; } - { - bool b_replicas = false; - get_keyval (conf, "multipleReplicas", b_replicas, false); - if (b_replicas) - comm = multiple_replicas; - else - comm = single_replica; - } - if (comm != single_replica) { if (expand_grids) cvm::fatal_error ("Error: expandBoundaries is not supported when " "using more than one replicas; please allocate " "wide enough boundaries for each colvar" "ahead of time.\n"); if (get_keyval (conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) { if (dump_replica_fes && (! dump_fes)) { cvm::log ("Enabling \"dumpFreeEnergyFile\".\n"); } } get_keyval (conf, "replicaID", replica_id, std::string ("")); if (!replica_id.size()) cvm::fatal_error ("Error: replicaID must be defined " "when using more than one replica.\n"); get_keyval (conf, "replicasRegistry", replicas_registry_file, (this->name+".replicas.registry.txt")); get_keyval (conf, "replicaUpdateFrequency", replica_update_freq, new_hill_freq); if (keep_hills) cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": keepHills with more than one replica can lead to a very " "large amount input/output and slow down your calculations. " "Please consider disabling it.\n"); { // TODO: one may want to specify the path manually for intricated filesystems? char *pwd = new char[3001]; if (GETCWD (pwd, 3000) == NULL) cvm::fatal_error ("Error: cannot get the path of the current working directory.\n"); replica_list_file = (std::string (pwd)+std::string (PATHSEP)+ this->name+"."+replica_id+".files.txt"); // replica_hills_file and replica_state_file are those written // by the current replica; within the mirror biases, they are // those by another replica replica_hills_file = (std::string (pwd)+std::string (PATHSEP)+ cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills"); replica_state_file = (std::string (pwd)+std::string (PATHSEP)+ cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state"); delete pwd; } // now register this replica // first check that it isn't already there bool registered_replica = false; std::ifstream reg_is (replicas_registry_file.c_str()); if (reg_is.good()) { // the file may not be there yet std::string existing_replica (""); std::string existing_replica_file (""); while ((reg_is >> existing_replica) && existing_replica.size() && (reg_is >> existing_replica_file) && existing_replica_file.size()) { if (existing_replica == replica_id) { // this replica was already registered replica_list_file = existing_replica_file; reg_is.close(); registered_replica = true; break; } } reg_is.close(); } // if this replica was not included yet, we should generate a // new record for it: but first, we write this replica's files, // for the others to read // open the "hills" buffer file replica_hills_os.open (replica_hills_file.c_str()); if (!replica_hills_os.good()) cvm::fatal_error ("Error: in opening file \""+ replica_hills_file+"\" for writing.\n"); replica_hills_os.setf (std::ios::scientific, std::ios::floatfield); // write the state file (so that there is always one available) write_replica_state_file(); // schedule to read the state files of the other replicas for (size_t ir = 0; ir < replicas.size(); ir++) { (replicas[ir])->replica_state_file_in_sync = false; } // if we're running without grids, use a growing list of "hills" files // otherwise, just one state file and one "hills" file as buffer std::ofstream list_os (replica_list_file.c_str(), (use_grids ? std::ios::trunc : std::ios::app)); if (! list_os.good()) cvm::fatal_error ("Error: in opening file \""+ replica_list_file+"\" for writing.\n"); list_os << "stateFile " << replica_state_file << "\n"; list_os << "hillsFile " << replica_hills_file << "\n"; list_os.close(); // finally, if add a new record for this replica to the registry if (! registered_replica) { std::ofstream reg_os (replicas_registry_file.c_str(), std::ios::app); if (! reg_os.good()) cvm::fatal_error ("Error: in opening file \""+ replicas_registry_file+"\" for writing.\n"); reg_os << replica_id << " " << replica_list_file << "\n"; reg_os.close(); } } get_keyval (conf, "writeHillsTrajectory", b_hills_traj, false); if (b_hills_traj) { std::string const traj_file_name (cvm::output_prefix+ ".colvars."+this->name+ ( (comm != single_replica) ? ("."+replica_id) : ("") )+ ".hills.traj"); hills_traj_os.open (traj_file_name.c_str()); if (!hills_traj_os.good()) cvm::fatal_error ("Error: in opening hills output file \"" + traj_file_name + "\".\n"); } // for well-tempered metadynamics get_keyval (conf, "wellTempered", well_tempered, false); get_keyval (conf, "biasTemperature", bias_temperature, -1.0); if ((bias_temperature == -1.0) && well_tempered) { cvm::fatal_error ("Error: biasTemperature is not set.\n"); } if (well_tempered) { cvm::log("Well-tempered metadynamics is used.\n"); cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n"); } if (cvm::debug()) cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); save_delimiters = false; } colvarbias_meta::~colvarbias_meta() { if (hills_energy) { delete hills_energy; hills_energy = NULL; } if (hills_energy_gradients) { delete hills_energy_gradients; hills_energy_gradients = NULL; } if (replica_hills_os.good()) replica_hills_os.close(); if (hills_traj_os.good()) hills_traj_os.close(); } // ********************************************************************** // Hill management member functions // ********************************************************************** std::list::const_iterator colvarbias_meta::create_hill (colvarbias_meta::hill const &h) { hill_iter const hills_end = hills.end(); hills.push_back (h); if (new_hills_begin == hills_end) { // if new_hills_begin is unset, set it for the first time new_hills_begin = hills.end(); new_hills_begin--; } if (use_grids) { // also add it to the list of hills that are off-grid, which may // need to be computed analytically when the colvar returns // off-grid cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers, true); if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) { hills_off_grid.push_back (h); } } // output to trajectory (if specified) if (hills_traj_os.good()) { hills_traj_os << (hills.back()).output_traj(); if (cvm::debug()) { hills_traj_os.flush(); } } has_data = true; return hills.end(); } std::list::const_iterator colvarbias_meta::delete_hill (hill_iter &h) { if (cvm::debug()) { cvm::log ("Deleting hill from the metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ", with step number "+ cvm::to_str (h->it)+(h->replica.size() ? ", replica id \""+h->replica : "")+".\n"); } if (use_grids && hills_off_grid.size()) { for (hill_iter hoff = hills_off_grid.begin(); hoff != hills_off_grid.end(); hoff++) { if (*h == *hoff) { hills_off_grid.erase (hoff); break; } } } if (hills_traj_os.good()) { // output to the trajectory hills_traj_os << "# DELETED this hill: " << (hills.back()).output_traj() << "\n"; if (cvm::debug()) hills_traj_os.flush(); } return hills.erase (h); } cvm::real colvarbias_meta::update() { if (cvm::debug()) cvm::log ("Updating the metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); if (use_grids) { std::vector curr_bin = hills_energy->get_colvars_index(); if (expand_grids) { // first of all, expand the grids, if specified if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": current coordinates on the grid: "+ cvm::to_str (curr_bin)+".\n"); bool changed_grids = false; int const min_buffer = (3 * (size_t) std::floor (hill_width)) + 1; std::vector new_sizes (hills_energy->sizes()); std::vector new_lower_boundaries (hills_energy->lower_boundaries); std::vector new_upper_boundaries (hills_energy->upper_boundaries); for (size_t i = 0; i < colvars.size(); i++) { if (! colvars[i]->expand_boundaries) continue; cvm::real &new_lb = new_lower_boundaries[i].real_value; cvm::real &new_ub = new_upper_boundaries[i].real_value; int &new_size = new_sizes[i]; bool changed_lb = false, changed_ub = false; - if (curr_bin[i] < min_buffer) { - int const extra_points = (min_buffer - curr_bin[i]); - new_lb -= extra_points * colvars[i]->width; - new_size += extra_points; - // changed offset in this direction => the pointer needs to - // be changed, too - curr_bin[i] += extra_points; + if (!colvars[i]->hard_lower_boundary) + if (curr_bin[i] < min_buffer) { + int const extra_points = (min_buffer - curr_bin[i]); + new_lb -= extra_points * colvars[i]->width; + new_size += extra_points; + // changed offset in this direction => the pointer needs to + // be changed, too + curr_bin[i] += extra_points; - changed_lb = true; - cvm::log ("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": new lower boundary for colvar \""+ - colvars[i]->name+"\", at "+ - cvm::to_str (new_lower_boundaries[i])+".\n"); - } + changed_lb = true; + cvm::log ("Metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": new lower boundary for colvar \""+ + colvars[i]->name+"\", at "+ + cvm::to_str (new_lower_boundaries[i])+".\n"); + } - if (curr_bin[i] > new_size - min_buffer - 1) { - int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer); - new_ub += extra_points * colvars[i]->width; - new_size += extra_points; + if (!colvars[i]->hard_upper_boundary) + if (curr_bin[i] > new_size - min_buffer - 1) { + int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer); + new_ub += extra_points * colvars[i]->width; + new_size += extra_points; - changed_ub = true; - cvm::log ("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": new upper boundary for colvar \""+ - colvars[i]->name+"\", at "+ - cvm::to_str (new_upper_boundaries[i])+".\n"); - } + changed_ub = true; + cvm::log ("Metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": new upper boundary for colvar \""+ + colvars[i]->name+"\", at "+ + cvm::to_str (new_upper_boundaries[i])+".\n"); + } if (changed_lb || changed_ub) changed_grids = true; } if (changed_grids) { // map everything into new grids colvar_grid_scalar *new_hills_energy = new colvar_grid_scalar (*hills_energy); colvar_grid_gradient *new_hills_energy_gradients = new colvar_grid_gradient (*hills_energy_gradients); // supply new boundaries to the new grids new_hills_energy->lower_boundaries = new_lower_boundaries; new_hills_energy->upper_boundaries = new_upper_boundaries; new_hills_energy->create (new_sizes, 0.0, 1); new_hills_energy_gradients->lower_boundaries = new_lower_boundaries; new_hills_energy_gradients->upper_boundaries = new_upper_boundaries; new_hills_energy_gradients->create (new_sizes, 0.0, colvars.size()); new_hills_energy->map_grid (*hills_energy); new_hills_energy_gradients->map_grid (*hills_energy_gradients); delete hills_energy; delete hills_energy_gradients; hills_energy = new_hills_energy; hills_energy_gradients = new_hills_energy_gradients; curr_bin = hills_energy->get_colvars_index(); if (cvm::debug()) cvm::log ("Coordinates on the new grid: "+ cvm::to_str (curr_bin)+".\n"); } } } // add a new hill if the required time interval has passed if ((cvm::step_absolute() % new_hill_freq) == 0) { if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": adding a new hill at step "+cvm::to_str (cvm::step_absolute())+".\n"); switch (comm) { case single_replica: if (well_tempered) { std::vector curr_bin = hills_energy->get_colvars_index(); cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width)); } else { create_hill (hill (hill_weight, colvars, hill_width)); } break; case multiple_replicas: if (well_tempered) { std::vector curr_bin = hills_energy->get_colvars_index(); cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id)); } else { create_hill (hill (hill_weight, colvars, hill_width, replica_id)); } if (replica_hills_os.good()) { replica_hills_os << hills.back(); } else { cvm::fatal_error ("Error: in metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ " while writing hills for the other replicas.\n"); } break; } } // sync with the other replicas (if needed) if (comm != single_replica) { // reread the replicas registry if ((cvm::step_absolute() % replica_update_freq) == 0) { update_replicas_registry(); // empty the output buffer replica_hills_os.flush(); read_replica_files(); } } // calculate the biasing energy and forces bias_energy = 0.0; for (size_t ir = 0; ir < colvars.size(); ir++) { colvar_forces[ir].reset(); } if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { replicas[ir]->bias_energy = 0.0; for (size_t ic = 0; ic < colvars.size(); ic++) { replicas[ir]->colvar_forces[ic].reset(); } } if (use_grids) { // get the forces from the grid if ((cvm::step_absolute() % grids_freq) == 0) { // map the most recent gaussians to the grids project_hills (new_hills_begin, hills.end(), hills_energy, hills_energy_gradients); new_hills_begin = hills.end(); // TODO: we may want to condense all into one replicas array, // including "this" as the first element if (comm == multiple_replicas) { for (size_t ir = 0; ir < replicas.size(); ir++) { replicas[ir]->project_hills (replicas[ir]->new_hills_begin, replicas[ir]->hills.end(), replicas[ir]->hills_energy, replicas[ir]->hills_energy_gradients); replicas[ir]->new_hills_begin = replicas[ir]->hills.end(); } } } std::vector curr_bin = hills_energy->get_colvars_index(); if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": current coordinates on the grid: "+ cvm::to_str (curr_bin)+".\n"); if (hills_energy->index_ok (curr_bin)) { // within the grid: add the energy and the forces from there bias_energy += hills_energy->value (curr_bin); for (size_t ic = 0; ic < colvars.size(); ic++) { cvm::real const *f = &(hills_energy_gradients->value (curr_bin)); colvar_forces[ic].real_value += -1.0 * f[ic]; // the gradients are stored, not the forces } if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { bias_energy += replicas[ir]->hills_energy->value (curr_bin); cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value (curr_bin)); for (size_t ic = 0; ic < colvars.size(); ic++) { colvar_forces[ic].real_value += -1.0 * f[ic]; } } } else { // off the grid: compute analytically only the hills at the grid's edges calc_hills (hills_off_grid.begin(), hills_off_grid.end(), bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force (ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces); } if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { calc_hills (replicas[ir]->hills_off_grid.begin(), replicas[ir]->hills_off_grid.end(), bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force (ic, replicas[ir]->hills_off_grid.begin(), replicas[ir]->hills_off_grid.end(), colvar_forces); } } } } // now include the hills that have not been binned yet (starting // from new_hills_begin) calc_hills (new_hills_begin, hills.end(), bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces); } if (cvm::debug()) cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+ ", hills forces = "+cvm::to_str (colvar_forces)+".\n"); if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": adding the forces from the other replicas.\n"); if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { calc_hills (replicas[ir]->new_hills_begin, replicas[ir]->hills.end(), bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force (ic, replicas[ir]->new_hills_begin, replicas[ir]->hills.end(), colvar_forces); } if (cvm::debug()) cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+ ", hills forces = "+cvm::to_str (colvar_forces)+".\n"); } return bias_energy; } void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, cvm::real &energy, std::vector const &colvar_values) { std::vector curr_values (colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { curr_values[i].type (colvars[i]->type()); } if (colvar_values.size()) { for (size_t i = 0; i < colvars.size(); i++) { curr_values[i] = colvar_values[i]; } } else { for (size_t i = 0; i < colvars.size(); i++) { curr_values[i] = colvars[i]->value(); } } for (hill_iter h = h_first; h != h_last; h++) { // compute the gaussian exponent cvm::real cv_sqdev = 0.0; for (size_t i = 0; i < colvars.size(); i++) { colvarvalue const &x = curr_values[i]; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; cv_sqdev += (colvars[i]->dist2 (x, center)) / (half_width*half_width); } // compute the gaussian if (cv_sqdev > 23.0) { // set it to zero if the exponent is more negative than log(1.0E-05) h->value (0.0); } else { h->value (std::exp (-0.5*cv_sqdev)); } energy += h->energy(); } } void colvarbias_meta::calc_hills_force (size_t const &i, colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, std::vector &forces, std::vector const &values) { // Retrieve the value of the colvar colvarvalue x (values.size() ? values[i].type() : colvars[i]->type()); x = (values.size() ? values[i] : colvars[i]->value()); // do the type check only once (all colvarvalues in the hills series // were already saved with their types matching those in the // colvars) switch (colvars[i]->type()) { case colvarvalue::type_scalar: for (hill_iter h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; forces[i].real_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * (colvars[i]->dist2_lgrad (x, center)).real_value ); } break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: for (hill_iter h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; forces[i].rvector_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * (colvars[i]->dist2_lgrad (x, center)).rvector_value ); } break; case colvarvalue::type_quaternion: for (hill_iter h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; forces[i].quaternion_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * (colvars[i]->dist2_lgrad (x, center)).quaternion_value ); } break; case colvarvalue::type_notset: break; } } // ********************************************************************** // grid management functions // ********************************************************************** void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, colvar_grid_scalar *he, colvar_grid_gradient *hg, bool print_progress) { if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": projecting hills.\n"); // TODO: improve it by looping over a small subgrid instead of the whole grid std::vector colvar_values (colvars.size()); std::vector colvar_forces_scalar (colvars.size()); std::vector he_ix = he->new_index(); std::vector hg_ix = (hg != NULL) ? hg->new_index() : std::vector (0); cvm::real hills_energy_here = 0.0; std::vector hills_forces_here (colvars.size(), 0.0); size_t count = 0; size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1))); if (hg != NULL) { // loop over the points of the grid for ( ; (he->index_ok (he_ix)) && (hg->index_ok (hg_ix)); count++) { for (size_t i = 0; i < colvars.size(); i++) { colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i); } // loop over the hills and increment the energy grid locally hills_energy_here = 0.0; calc_hills (h_first, h_last, hills_energy_here, colvar_values); he->acc_value (he_ix, hills_energy_here); for (size_t i = 0; i < colvars.size(); i++) { hills_forces_here[i].reset(); calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values); colvar_forces_scalar[i] = hills_forces_here[i].real_value; } hg->acc_force (hg_ix, &(colvar_forces_scalar.front())); he->incr (he_ix); hg->incr (hg_ix); if ((count % print_frequency) == 0) { if (print_progress) { cvm::real const progress = cvm::real (count) / cvm::real (hg->number_of_points()); std::ostringstream os; os.setf (std::ios::fixed, std::ios::floatfield); os << std::setw (6) << std::setprecision (2) << 100.0 * progress << "% done."; cvm::log (os.str()); } } } } else { // simpler version, with just the energy for ( ; (he->index_ok (he_ix)); ) { for (size_t i = 0; i < colvars.size(); i++) { colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i); } hills_energy_here = 0.0; calc_hills (h_first, h_last, hills_energy_here, colvar_values); he->acc_value (he_ix, hills_energy_here); he->incr (he_ix); count++; if ((count % print_frequency) == 0) { if (print_progress) { cvm::real const progress = cvm::real (count) / cvm::real (he->number_of_points()); std::ostringstream os; os.setf (std::ios::fixed, std::ios::floatfield); os << std::setw (6) << std::setprecision (2) << 100.0 * progress << "% done."; cvm::log (os.str()); } } } } if (print_progress) { cvm::log ("100.00% done."); } if (! keep_hills) { hills.erase (hills.begin(), hills.end()); } } void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, colvar_grid_scalar *he) { hills_off_grid.clear(); for (hill_iter h = h_first; h != h_last; h++) { cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers, true); if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) { hills_off_grid.push_back (*h); } } } // ********************************************************************** // multiple replicas functions // ********************************************************************** void colvarbias_meta::update_replicas_registry() { if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": updating the list of replicas, currently containing "+ cvm::to_str (replicas.size())+" elements.\n"); { // copy the whole file into a string for convenience std::string line (""); std::ifstream reg_file (replicas_registry_file.c_str()); if (reg_file.good()) { replicas_registry.clear(); while (colvarparse::getline_nocomments (reg_file, line)) replicas_registry.append (line+"\n"); } else { cvm::fatal_error ("Error: failed to open file \""+replicas_registry_file+ "\" for reading.\n"); } } // now parse it std::istringstream reg_is (replicas_registry); if (reg_is.good()) { std::string new_replica (""); std::string new_replica_file (""); while ((reg_is >> new_replica) && new_replica.size() && (reg_is >> new_replica_file) && new_replica_file.size()) { if (new_replica == this->replica_id) { // this is the record for this same replica, skip it if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": skipping this replica's own record: \""+ new_replica+"\", \""+new_replica_file+"\"\n"); new_replica_file.clear(); new_replica.clear(); continue; } bool already_loaded = false; for (size_t ir = 0; ir < replicas.size(); ir++) { if (new_replica == (replicas[ir])->replica_id) { // this replica was already added if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": skipping a replica already loaded, \""+ (replicas[ir])->replica_id+"\".\n"); already_loaded = true; break; } } if (!already_loaded) { // add this replica to the registry cvm::log ("Metadynamics bias \""+this->name+"\""+ ": accessing replica \""+new_replica+"\".\n"); replicas.push_back (new colvarbias_meta()); (replicas.back())->replica_id = new_replica; (replicas.back())->replica_list_file = new_replica_file; (replicas.back())->replica_state_file = ""; (replicas.back())->replica_state_file_in_sync = false; // Note: the following could become a copy constructor? (replicas.back())->name = this->name; (replicas.back())->colvars = colvars; (replicas.back())->use_grids = use_grids; (replicas.back())->dump_fes = false; (replicas.back())->expand_grids = false; (replicas.back())->rebin_grids = false; (replicas.back())->keep_hills = false; (replicas.back())->colvar_forces = colvar_forces; (replicas.back())->comm = multiple_replicas; if (use_grids) { (replicas.back())->hills_energy = new colvar_grid_scalar (colvars); (replicas.back())->hills_energy_gradients = new colvar_grid_gradient (colvars); } } } // continue the cycle new_replica_file = ""; new_replica = ""; } else { cvm::fatal_error ("Error: cannot read the replicas registry file \""+ replicas_registry+"\".\n"); } // now (re)read the list file of each replica for (size_t ir = 0; ir < replicas.size(); ir++) { if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": reading the list file for replica \""+ (replicas[ir])->replica_id+"\".\n"); std::ifstream list_is ((replicas[ir])->replica_list_file.c_str()); std::string key; std::string new_state_file, new_hills_file; if (!(list_is >> key) || !(key == std::string ("stateFile")) || !(list_is >> new_state_file) || !(list_is >> key) || !(key == std::string ("hillsFile")) || !(list_is >> new_hills_file)) { cvm::log ("Metadynamics bias \""+this->name+"\""+ ": failed to read the file \""+ (replicas[ir])->replica_list_file+"\": will try again after "+ cvm::to_str (replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; } else { (replicas[ir])->update_status = 0; if (new_state_file != (replicas[ir])->replica_state_file) { cvm::log ("Metadynamics bias \""+this->name+"\""+ ": replica \""+(replicas[ir])->replica_id+ "\" has supplied a new state file, \""+new_state_file+ "\".\n"); (replicas[ir])->replica_state_file_in_sync = false; (replicas[ir])->replica_state_file = new_state_file; (replicas[ir])->replica_hills_file = new_hills_file; } } } if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\": the list of replicas contains "+ cvm::to_str (replicas.size())+" elements.\n"); } void colvarbias_meta::read_replica_files() { for (size_t ir = 0; ir < replicas.size(); ir++) { if (! (replicas[ir])->replica_state_file_in_sync) { // if a new state file is being read, the hills file is also new (replicas[ir])->replica_hills_file_pos = 0; } // (re)read the state file if necessary if ( (! (replicas[ir])->has_data) || (! (replicas[ir])->replica_state_file_in_sync) ) { cvm::log ("Metadynamics bias \""+this->name+"\""+ ": reading the state of replica \""+ (replicas[ir])->replica_id+"\" from file \""+ (replicas[ir])->replica_state_file+"\".\n"); std::ifstream is ((replicas[ir])->replica_state_file.c_str()); if (! (replicas[ir])->read_restart (is)) { cvm::log ("Reading from file \""+(replicas[ir])->replica_state_file+ "\" failed or incomplete: will try again in "+ cvm::to_str (replica_update_freq)+" steps.\n"); } else { // state file has been read successfully (replicas[ir])->replica_state_file_in_sync = true; (replicas[ir])->update_status = 0; } is.close(); } // now read the hills added after writing the state file if ((replicas[ir])->replica_hills_file.size()) { if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": checking for new hills from replica \""+ (replicas[ir])->replica_id+"\" in the file \""+ (replicas[ir])->replica_hills_file+"\".\n"); // read hills from the other replicas' files; for each file, resume // the position recorded previously std::ifstream is ((replicas[ir])->replica_hills_file.c_str()); if (is.good()) { // try to resume the previous position is.seekg ((replicas[ir])->replica_hills_file_pos, std::ios::beg); if (!is.good()){ // if fail (the file may have been overwritten), reset this // position is.clear(); is.seekg (0, std::ios::beg); // reset the counter (replicas[ir])->replica_hills_file_pos = 0; // schedule to reread the state file (replicas[ir])->replica_state_file_in_sync = false; // and record the failure (replicas[ir])->update_status++; cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+ "\" at the previous position: will try again in "+ cvm::to_str (replica_update_freq)+" steps.\n"); } else { while ((replicas[ir])->read_hill (is)) { // if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": received a hill from replica \""+ (replicas[ir])->replica_id+ "\" at step "+ cvm::to_str (((replicas[ir])->hills.back()).it)+".\n"); } is.clear(); // store the position for the next read (replicas[ir])->replica_hills_file_pos = is.tellg(); if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": stopped reading file \""+(replicas[ir])->replica_hills_file+ "\" at position "+ cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n"); // test whether this is the end of the file is.seekg (0, std::ios::end); if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) { (replicas[ir])->update_status++; } else { (replicas[ir])->update_status = 0; } } } else { cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+ "\": will try again in "+ cvm::to_str (replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; // cvm::fatal_error ("Error: cannot read from file \""+ // (replicas[ir])->replica_hills_file+"\".\n"); } is.close(); } size_t const n_flush = (replica_update_freq/new_hill_freq + 1); if ((replicas[ir])->update_status > 3*n_flush) { // TODO: suspend the calculation? cvm::log ("WARNING: in metadynamics bias \""+this->name+"\""+ " failed to read completely the output of replica \""+ (replicas[ir])->replica_id+ "\" after more than "+ cvm::to_str ((replicas[ir])->update_status * replica_update_freq)+ " steps. Ensure that it is still running.\n"); } } } // ********************************************************************** // input functions // ********************************************************************** std::istream & colvarbias_meta::read_restart (std::istream& is) { size_t const start_pos = is.tellg(); if (comm == single_replica) { // if using a multiple replicas scheme, output messages // are printed before and after calling this function cvm::log ("Restarting metadynamics bias \""+this->name+"\""+ ".\n"); } std::string key, brace, conf; if ( !(is >> key) || !(key == "metadynamics") || !(is >> brace) || !(brace == "{") || !(is >> colvarparse::read_block ("configuration", conf)) ) { if (comm == single_replica) cvm::log ("Error: in reading restart configuration for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ (replica_state_file_in_sync ? ("at position "+ cvm::to_str (start_pos)+ " in the state file") : "")+".\n"); is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } std::string name = ""; if ( colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent) && (name != this->name) ) cvm::fatal_error ("Error: in the restart file, the " "\"metadynamics\" block has a different name: different system?\n"); if (name.size() == 0) { cvm::fatal_error ("Error: \"metadynamics\" block within the restart file " "has no identifiers.\n"); } if (comm != single_replica) { std::string replica = ""; if (colvarparse::get_keyval (conf, "replicaID", replica, std::string (""), colvarparse::parse_silent) && (replica != this->replica_id)) cvm::fatal_error ("Error: in the restart file, the " "\"metadynamics\" block has a different replicaID: different system?\n"); colvarparse::get_keyval (conf, "step", state_file_step, cvm::step_absolute(), colvarparse::parse_silent); } bool grids_from_restart_file = use_grids; if (use_grids) { if (expand_grids) { // the boundaries of the colvars may have been changed; TODO: // this reallocation is only for backward-compatibility, and may // be deleted when grid_parameters (i.e. colvargrid's own // internal reallocation) has kicked in delete hills_energy; delete hills_energy_gradients; hills_energy = new colvar_grid_scalar (colvars); hills_energy_gradients = new colvar_grid_gradient (colvars); } colvar_grid_scalar *hills_energy_backup = NULL; colvar_grid_gradient *hills_energy_gradients_backup = NULL; if (has_data) { if (cvm::debug()) cvm::log ("Backupping grids for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); hills_energy_backup = hills_energy; hills_energy_gradients_backup = hills_energy_gradients; hills_energy = new colvar_grid_scalar (colvars); hills_energy_gradients = new colvar_grid_gradient (colvars); } size_t const hills_energy_pos = is.tellg(); if (!(is >> key)) { if (hills_energy_backup != NULL) { delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; hills_energy_gradients = hills_energy_gradients_backup; } is.clear(); is.seekg (hills_energy_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } else if (!(key == std::string ("hills_energy")) || !(hills_energy->read_restart (is))) { is.clear(); is.seekg (hills_energy_pos, std::ios::beg); grids_from_restart_file = false; if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error ("Error: couldn't read the free energy grid for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ "; if useGrids was off when the state file was written, " "enable rebinGrids now to regenerate the grids.\n"); else { if (comm == single_replica) cvm::log ("Error: couldn't read the free energy grid for metadynamics bias \""+ this->name+"\".\n"); delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; hills_energy_gradients = hills_energy_gradients_backup; is.setstate (std::ios::failbit); return is; } } } size_t const hills_energy_gradients_pos = is.tellg(); if (!(is >> key)) { if (hills_energy_backup != NULL) { delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; hills_energy_gradients = hills_energy_gradients_backup; } is.clear(); is.seekg (hills_energy_gradients_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } else if (!(key == std::string ("hills_energy_gradients")) || !(hills_energy_gradients->read_restart (is))) { is.clear(); is.seekg (hills_energy_gradients_pos, std::ios::beg); grids_from_restart_file = false; if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ "; if useGrids was off when the state file was written, " "enable rebinGrids now to regenerate the grids.\n"); else { if (comm == single_replica) cvm::log ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+ this->name+"\".\n"); delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; hills_energy_gradients = hills_energy_gradients_backup; is.setstate (std::ios::failbit); return is; } } } if (cvm::debug()) cvm::log ("Successfully read new grids for bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); if (hills_energy_backup != NULL) { // now that we have successfully updated the grids, delete the // backup copies if (cvm::debug()) cvm::log ("Deallocating the older grids.\n"); delete hills_energy_backup; delete hills_energy_gradients_backup; } } bool const existing_hills = (hills.size() > 0); size_t const old_hills_size = hills.size(); hill_iter old_hills_end = hills.end(); hill_iter old_hills_off_grid_end = hills_off_grid.end(); // read the hills explicitly written (if there are any) while (read_hill (is)) { if (cvm::debug()) cvm::log ("Read a previously saved hill under the " "metadynamics bias \""+ this->name+"\", created at step "+ cvm::to_str ((hills.back()).it)+".\n"); } is.clear(); new_hills_begin = hills.end(); if (grids_from_restart_file) { if (hills.size() > old_hills_size) cvm::log ("Read "+cvm::to_str (hills.size())+ " hills in addition to the grids.\n"); } else { if (hills.size()) cvm::log ("Read "+cvm::to_str (hills.size())+" hills.\n"); } if (rebin_grids) { // allocate new grids (based on the new boundaries and widths just // read from the configuration file), and project onto them the // grids just read from the restart file colvar_grid_scalar *new_hills_energy = new colvar_grid_scalar (colvars); colvar_grid_gradient *new_hills_energy_gradients = new colvar_grid_gradient (colvars); if (!grids_from_restart_file || (keep_hills && hills.size())) { // if there are hills, recompute the new grids from them cvm::log ("Rebinning the energy and forces grids from "+ cvm::to_str (hills.size())+" hills (this may take a while)...\n"); project_hills (hills.begin(), hills.end(), new_hills_energy, new_hills_energy_gradients, true); cvm::log ("rebinning done.\n"); } else { // otherwise, use the grids in the restart file cvm::log ("Rebinning the energy and forces grids " "from the grids in the restart file.\n"); new_hills_energy->map_grid (*hills_energy); new_hills_energy_gradients->map_grid (*hills_energy_gradients); } delete hills_energy; delete hills_energy_gradients; hills_energy = new_hills_energy; hills_energy_gradients = new_hills_energy_gradients; // assuming that some boundaries have expanded, eliminate those // off-grid hills that aren't necessary any more if (hills.size()) recount_hills_off_grid (hills.begin(), hills.end(), hills_energy); } if (use_grids) { if (hills_off_grid.size()) { cvm::log (cvm::to_str (hills_off_grid.size())+" hills are near the " "grid boundaries: they will be computed analytically " "and saved to the state files.\n"); } } is >> brace; if (brace != "}") { cvm::log ("Incomplete restart information for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": no closing brace at position "+ cvm::to_str (is.tellg())+" in the file.\n"); is.setstate (std::ios::failbit); return is; } if (cvm::debug()) cvm::log ("colvarbias_meta::read_restart() done\n"); if (existing_hills) { hills.erase (hills.begin(), old_hills_end); hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end); } has_data = true; if (comm != single_replica) { read_replica_files(); } return is; } std::istream & colvarbias_meta::read_hill (std::istream &is) { if (!is) return is; // do nothing if failbit is set size_t const start_pos = is.tellg(); std::string data; if ( !(is >> read_block ("hill", data)) ) { is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } size_t h_it; get_keyval (data, "step", h_it, 0, parse_silent); if (h_it <= state_file_step) { if (cvm::debug()) cvm::log ("Skipping a hill older than the state file for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); return is; } cvm::real h_weight; get_keyval (data, "weight", h_weight, hill_weight, parse_silent); std::vector h_centers (colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { h_centers[i].type ((colvars[i]->value()).type()); } { // it is safer to read colvarvalue objects one at a time; // TODO: change this it later std::string centers_input; key_lookup (data, "centers", centers_input); std::istringstream centers_is (centers_input); for (size_t i = 0; i < colvars.size(); i++) { centers_is >> h_centers[i]; } } std::vector h_widths (colvars.size()); get_keyval (data, "widths", h_widths, std::vector (colvars.size(), (std::sqrt (2.0 * PI) / 2.0)), parse_silent); std::string h_replica = ""; if (comm != single_replica) { get_keyval (data, "replicaID", h_replica, replica_id, parse_silent); if (h_replica != replica_id) cvm::fatal_error ("Error: trying to read a hill created by replica \""+h_replica+ "\" for replica \""+replica_id+ "\"; did you swap output files?\n"); } hill_iter const hills_end = hills.end(); hills.push_back (hill (h_it, h_weight, h_centers, h_widths, h_replica)); if (new_hills_begin == hills_end) { // if new_hills_begin is unset, set it for the first time new_hills_begin = hills.end(); new_hills_begin--; } if (use_grids) { // add this also to the list of hills that are off-grid, which will // be computed analytically cvm::real const min_dist = hills_energy->bin_distance_from_boundaries ((hills.back()).centers, true); if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) { hills_off_grid.push_back (hills.back()); } } has_data = true; return is; } // ********************************************************************** // output functions // ********************************************************************** std::ostream & colvarbias_meta::write_restart (std::ostream& os) { os << "metadynamics {\n" << " configuration {\n" << " step " << cvm::step_absolute() << "\n" << " name " << this->name << "\n"; if (this->comm != single_replica) os << " replicaID " << this->replica_id << "\n"; os << " }\n\n"; if (use_grids) { // this is a very good time to project hills, if you haven't done // it already! project_hills (new_hills_begin, hills.end(), hills_energy, hills_energy_gradients); new_hills_begin = hills.end(); // write down the grids to the restart file os << " hills_energy\n"; hills_energy->write_restart (os); os << " hills_energy_gradients\n"; hills_energy_gradients->write_restart (os); } if ( (!use_grids) || keep_hills ) { // write all hills currently in memory for (std::list::const_iterator h = this->hills.begin(); h != this->hills.end(); h++) { os << *h; } } else { // write just those that are near the grid boundaries for (std::list::const_iterator h = this->hills_off_grid.begin(); h != this->hills_off_grid.end(); h++) { os << *h; } } os << "}\n\n"; if (comm != single_replica) { write_replica_state_file(); // schedule to reread the state files of the other replicas (they // have also rewritten them) for (size_t ir = 0; ir < replicas.size(); ir++) { (replicas[ir])->replica_state_file_in_sync = false; } } if (dump_fes) { write_pmf(); } return os; } void colvarbias_meta::write_pmf() { // allocate a new grid to store the pmf - colvar_grid_scalar *pmf = new colvar_grid_scalar (colvars); + colvar_grid_scalar *pmf = new colvar_grid_scalar (*hills_energy); + pmf->create(); std::string fes_file_name_prefix (cvm::output_prefix); if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) { // if this is not the only free energy integrator, append // this bias's name, to distinguish it from the output of the other // biases producing a .pmf file // TODO: fix for ABF with updateBias == no fes_file_name_prefix += ("."+this->name); } if ((comm == single_replica) || (dump_replica_fes)) { // output the PMF from this instance or replica pmf->reset(); pmf->add_grid (*hills_energy); cvm::real const max = pmf->maximum_value(); pmf->add_constant (-1.0 * max); pmf->multiply_constant (-1.0); if (well_tempered) { cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature; pmf->multiply_constant (well_temper_scale); } - std::string const fes_file_name (fes_file_name_prefix + - ((comm != single_replica) ? ".partial" : "") + - (dump_fes_save ? - "."+cvm::to_str (cvm::step_absolute()) : "") + - ".pmf"); - cvm::backup_file (fes_file_name.c_str()); - std::ofstream fes_os (fes_file_name.c_str()); - pmf->write_multicol (fes_os); - fes_os.close(); + { + std::string const fes_file_name (fes_file_name_prefix + + ((comm != single_replica) ? ".partial" : "") + + (dump_fes_save ? + "."+cvm::to_str (cvm::step_absolute()) : "") + + ".pmf"); + cvm::backup_file (fes_file_name.c_str()); + std::ofstream fes_os (fes_file_name.c_str()); + pmf->write_multicol (fes_os); + fes_os.close(); + } } - if (comm != single_replica) { // output the combined PMF from all replicas pmf->reset(); pmf->add_grid (*hills_energy); for (size_t ir = 0; ir < replicas.size(); ir++) { pmf->add_grid (*(replicas[ir]->hills_energy)); } cvm::real const max = pmf->maximum_value(); pmf->add_constant (-1.0 * max); pmf->multiply_constant (-1.0); if (well_tempered) { cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature; pmf->multiply_constant (well_temper_scale); } std::string const fes_file_name (fes_file_name_prefix + (dump_fes_save ? "."+cvm::to_str (cvm::step_absolute()) : "") + ".pmf"); cvm::backup_file (fes_file_name.c_str()); std::ofstream fes_os (fes_file_name.c_str()); pmf->write_multicol (fes_os); fes_os.close(); } delete pmf; } void colvarbias_meta::write_replica_state_file() { // write down also the restart for the other replicas: TODO: this // is duplicated code, that could be cleaned up later cvm::backup_file (replica_state_file.c_str()); std::ofstream rep_state_os (replica_state_file.c_str()); if (!rep_state_os.good()) cvm::fatal_error ("Error: in opening file \""+ replica_state_file+"\" for writing.\n"); rep_state_os.setf (std::ios::scientific, std::ios::floatfield); rep_state_os << "\n" << "metadynamics {\n" << " configuration {\n" << " name " << this->name << "\n" << " step " << cvm::step_absolute() << "\n"; if (this->comm != single_replica) { rep_state_os << " replicaID " << this->replica_id << "\n"; } rep_state_os << " }\n\n"; rep_state_os << " hills_energy\n"; rep_state_os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width); hills_energy->write_restart (rep_state_os); rep_state_os << " hills_energy_gradients\n"; rep_state_os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width); hills_energy_gradients->write_restart (rep_state_os); if ( (!use_grids) || keep_hills ) { // write all hills currently in memory for (std::list::const_iterator h = this->hills.begin(); h != this->hills.end(); h++) { rep_state_os << *h; } } else { // write just those that are near the grid boundaries for (std::list::const_iterator h = this->hills_off_grid.begin(); h != this->hills_off_grid.end(); h++) { rep_state_os << *h; } } rep_state_os << "}\n\n"; rep_state_os.close(); // reopen the hills file replica_hills_os.close(); replica_hills_os.open (replica_hills_file.c_str()); if (!replica_hills_os.good()) cvm::fatal_error ("Error: in opening file \""+ replica_hills_file+"\" for writing.\n"); replica_hills_os.setf (std::ios::scientific, std::ios::floatfield); } std::string colvarbias_meta::hill::output_traj() { std::ostringstream os; os.setf (std::ios::fixed, std::ios::floatfield); os << std::setw (cvm::it_width) << it << " "; os.setf (std::ios::scientific, std::ios::floatfield); os << " "; for (size_t i = 0; i < centers.size(); i++) { os << " "; os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << centers[i]; } os << " "; for (size_t i = 0; i < widths.size(); i++) { os << " "; os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << widths[i]; } os << " "; os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << W << "\n"; return os.str(); } std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) { os.setf (std::ios::scientific, std::ios::floatfield); os << "hill {\n"; os << " step " << std::setw (cvm::it_width) << h.it << "\n"; os << " weight " << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << h.W << "\n"; if (h.replica.size()) os << " replicaID " << h.replica << "\n"; os << " centers "; for (size_t i = 0; i < (h.centers).size(); i++) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << h.centers[i]; } os << "\n"; os << " widths "; for (size_t i = 0; i < (h.widths).size(); i++) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << h.widths[i]; } os << "\n"; os << "}\n"; return os; } diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index f1997a6dd..b3e89527c 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -1,73 +1,146 @@ #include "colvarmodule.h" #include "colvarvalue.h" #include "colvar.h" #include "colvarcomp.h" colvar::cvc::cvc() : sup_coeff (1.0), sup_np (1), b_periodic (false), b_debug_gradients (false), b_inverse_gradients (false), b_Jacobian_derivative (false) {} colvar::cvc::cvc (std::string const &conf) : sup_coeff (1.0), sup_np (1), b_periodic (false), b_debug_gradients (false), b_inverse_gradients (false), b_Jacobian_derivative (false) { if (cvm::debug()) cvm::log ("Initializing cvc base object.\n"); get_keyval (conf, "name", this->name, std::string (""), parse_silent); get_keyval (conf, "componentCoeff", sup_coeff, 1.0); get_keyval (conf, "componentExp", sup_np, 1); get_keyval (conf, "period", period, 0.0); get_keyval (conf, "wrapAround", wrap_center, 0.0); get_keyval (conf, "debugGradients", b_debug_gradients, false, parse_silent); if (cvm::debug()) cvm::log ("Done initializing cvc base object.\n"); } void colvar::cvc::parse_group (std::string const &conf, char const *group_key, cvm::atom_group &group, bool optional) { if (key_lookup (conf, group_key)) { group.parse (conf, group_key); } else { if (! optional) { cvm::fatal_error ("Error: definition for atom group \""+ std::string (group_key)+"\" not found.\n"); } } } colvar::cvc::~cvc() {} void colvar::cvc::calc_force_invgrads() { cvm::fatal_error ("Error: calculation of inverse gradients is not implemented " "for colvar components of type \""+function_type+"\".\n"); } void colvar::cvc::calc_Jacobian_derivative() { cvm::fatal_error ("Error: calculation of inverse gradients is not implemented " "for colvar components of type \""+function_type+"\".\n"); } + + +void colvar::cvc::debug_gradients (cvm::atom_group &group) +{ + // this function should work for any scalar variable: + // the only difference will be the name of the atom group (here, "group") + + // collect into a vector for convenience + std::vector gradients (group.size()); + for (size_t i = 0; i < group.size(); i++) { + gradients[i] = group[i].grad; + } + + cvm::rotation const rot_0 = group.rot; + cvm::rotation const rot_inv = group.rot.inverse(); + + cvm::real const x_0 = x.real_value; + + cvm::log ("gradients = "+cvm::to_str (gradients)+"\n"); + + // it only makes sense to debug the fit gradients + // when the fitting group is the same as this group + if (group.b_fit_gradients && (group.ref_pos_group == NULL)) { + group.calc_fit_gradients(); + if (group.b_rotate) { + // fit_gradients are in the original frame, we should print them in the rotated frame + for (size_t j = 0; j < group.fit_gradients.size(); j++) { + group.fit_gradients[j] = rot_0.rotate (group.fit_gradients[j]); + } + } + cvm::log ("fit_gradients = "+cvm::to_str (group.fit_gradients)+"\n"); + if (group.b_rotate) { + for (size_t j = 0; j < group.fit_gradients.size(); j++) { + group.fit_gradients[j] = rot_inv.rotate (group.fit_gradients[j]); + } + } + } + + for (size_t ia = 0; ia < group.size(); ia++) { + + // tests are best conducted in the unrotated (simulation) frame + cvm::rvector const atom_grad = group.b_rotate ? + rot_inv.rotate (group[ia].grad) : + group[ia].grad; + + for (size_t id = 0; id < 3; id++) { + // (re)read original positions + group.read_positions(); + // change one coordinate + group[ia].pos[id] += cvm::debug_gradients_step_size; + // (re)do the fit (if defined) + if (group.b_center || group.b_rotate) { + group.calc_apply_roto_translation(); + } + calc_value(); + cvm::real const x_1 = x.real_value; + cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n"); + cvm::log ("dx(actual) = "+cvm::to_str (x_1 - x_0, + 21, 14)+"\n"); + cvm::real const dx_pred = (group.b_fit_gradients && (group.ref_pos_group == NULL)) ? + (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : + (cvm::debug_gradients_step_size * atom_grad[id]); + cvm::log ("dx(interp) = "+cvm::to_str (dx_pred, + 21, 14)+"\n"); + cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+ + cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) / + std::fabs (x_1 - x_0), + 12, 5)+ + ".\n"); + + } + } +} diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 33bd73f64..a8a06814f 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -1,1433 +1,1434 @@ -#ifndef COLVARDEF_H -#define COLVARDEF_H +#ifndef COLVARCOMP_H +#define COLVARCOMP_H #include #include #include "colvarmodule.h" #include "colvar.h" #include "colvaratoms.h" /// \brief Colvar component (base class); most implementations of /// \link cvc \endlink utilize one or more \link /// colvarmodule::atom \endlink or \link colvarmodule::atom_group /// \endlink objects to access atoms. /// /// A \link cvc \endlink object (or an object of a /// cvc-derived class) specifies how to calculate a collective /// variable, its gradients and other related physical quantities /// which do not depend only on the numeric value (the \link colvar /// \endlink class already serves this purpose). /// /// No restriction is set to what kind of calculation a \link /// cvc \endlink object performs (usually calculate an /// analytical function of atomic coordinates). The only constraint /// is that the value calculated is implemented as a \link colvarvalue /// \endlink object. This serves to provide a unique way to calculate /// scalar and non-scalar collective variables, and specify if and how /// they can be combined together by the parent \link colvar \endlink /// object. /// /// If you wish to implement a new collective variable component, you /// should write your own class by inheriting directly from \link /// cvc \endlink, or one of its derived classes (for instance, /// \link distance \endlink is frequently used, because it provides /// useful data and function members for any colvar based on two /// atom groups). The steps are: \par /// 1. add the name of this class under colvar.h \par /// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par /// 3. declare the class in colvarcomp.h \par /// 4. implement the class in one of the files colvarcomp_*.C /// /// /// The cvm::atom and cvm::atom_group classes are available to /// transparently communicate with the simulation program. However, /// they are not strictly needed, as long as all the degrees of /// freedom associated to the cvc are properly evolved by a simple /// call to e.g. apply_force(). class colvar::cvc : public colvarparse { public: /// \brief The name of the object (helps to identify this /// cvc instance when debugging) std::string name; /// \brief Description of the type of collective variable /// /// Normally this string is set by the parent \link colvar \endlink /// object within its constructor, when all \link cvc \endlink /// objects are initialized; therefore the main "config string" /// constructor does not need to define it. If a \link cvc /// \endlink is initialized and/or a different constructor is used, /// this variable definition should be set within the constructor. std::string function_type; /// \brief Type of \link colvarvalue \endlink that this cvc /// provides colvarvalue::Type type() const; /// \brief Coefficient in the polynomial combination (default: 1.0) cvm::real sup_coeff; /// \brief Exponent in the polynomial combination (default: 1) int sup_np; /// \brief Period of this cvc value, (default: 0.0, non periodic) cvm::real period; /// \brief If the component is periodic, wrap around this value (default: 0.0) cvm::real wrap_center; bool b_periodic; /// \brief Constructor /// /// At least one constructor which reads a string should be defined /// for every class inheriting from cvc \param conf Contents /// of the configuration file pertaining to this \link cvc /// \endlink cvc (std::string const &conf); /// \brief Within the constructor, make a group parse its own /// options from the provided configuration string void parse_group (std::string const &conf, char const *group_key, cvm::atom_group &group, bool optional = false); /// \brief Default constructor (used when \link cvc \endlink /// objects are declared within other ones) cvc(); /// Destructor virtual ~cvc(); - /// \brief If true, calc_gradients() will calculate - /// finite-difference gradients alongside the analytical ones and - /// report their differences - bool b_debug_gradients; /// \brief If this flag is false (default), inverse gradients /// (derivatives of atom coordinates with respect to x) are /// unavailable; it should be set to true by the constructor of each /// derived object capable of calculating them bool b_inverse_gradients; /// \brief If this flag is false (default), the Jacobian derivative /// (divergence of the inverse gradients) is unavailable; it should /// be set to true by the constructor of each derived object capable /// of calculating it bool b_Jacobian_derivative; /// \brief Calculate the variable virtual void calc_value() = 0; /// \brief Calculate the atomic gradients, to be reused later in /// order to apply forces virtual void calc_gradients() = 0; + /// \brief If true, calc_gradients() will call debug_gradients() for every group needed + bool b_debug_gradients; + + /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component + virtual void debug_gradients (cvm::atom_group &group); + /// \brief Calculate the total force from the system using the /// inverse atomic gradients virtual void calc_force_invgrads(); /// \brief Calculate the divergence of the inverse atomic gradients virtual void calc_Jacobian_derivative(); /// \brief Return the previously calculated value virtual colvarvalue value() const; /// \brief Return the previously calculated system force virtual colvarvalue system_force() const; /// \brief Return the previously calculated divergence of the /// inverse atomic gradients virtual colvarvalue Jacobian_derivative() const; /// \brief Apply the collective variable force, by communicating the /// atomic forces to the simulation program (\b Note: the \link ft /// \endlink member is not altered by this function) /// /// Note: multiple calls to this function within the same simulation /// step will add the forces altogether \param cvforce The /// collective variable force, usually coming from the biases and /// eventually manipulated by the parent \link colvar \endlink /// object virtual void apply_force (colvarvalue const &cvforce) = 0; /// \brief Square distance between x1 and x2 (can be redefined to /// transparently implement constraints, symmetries and /// periodicities) /// /// colvar::cvc::dist2() and the related functions are /// declared as "const" functions, but not "static", because /// additional parameters defining the metrics (e.g. the /// periodicity) may be specific to each colvar::cvc object. /// /// If symmetries or periodicities are present, the /// colvar::cvc::dist2() should be redefined to return the /// "closest distance" value and colvar::cvc::dist2_lgrad(), /// colvar::cvc::dist2_rgrad() to return its gradients. /// /// If constraints are present (and not already implemented by any /// of the \link colvarvalue \endlink types), the /// colvar::cvc::dist2_lgrad() and /// colvar::cvc::dist2_rgrad() functions should be redefined /// to provide a gradient which is compatible with the constraint, /// i.e. already deprived of its component normal to the constraint /// hypersurface. /// /// Finally, another useful application, if you are performing very /// many operations with these functions, could be to override the /// \link colvarvalue \endlink member functions and access directly /// its member data. For instance: to define dist2(x1,x2) as /// (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in /// case of a scalar \link colvarvalue \endlink type. virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Gradient (with respect to x1) of the square distance (can /// be redefined to transparently implement constraints, symmetries /// and periodicities) virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Gradient (with respect to x2) of the square distance (can /// be redefined to transparently implement constraints, symmetries /// and periodicities) virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Return a positive number if x2>x1, zero if x2==x1, /// negative otherwise (can be redefined to transparently implement /// constraints, symmetries and periodicities) \b Note: \b it \b /// only \b works \b with \b scalar \b variables, otherwise raises /// an error virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Wrapp value (for periodic/symmetric cvcs) virtual void wrap (colvarvalue &x) const; /// \brief Pointers to all atom groups, to let colvars collect info /// e.g. atomic gradients std::vector atom_groups; protected: /// \brief Cached value colvarvalue x; /// \brief Value at the previous step colvarvalue x_old; /// \brief Calculated system force (\b Note: this is calculated from /// the total atomic forces read from the program, subtracting fromt /// the "internal" forces of the system the "external" forces from /// the colvar biases) colvarvalue ft; /// \brief Calculated Jacobian derivative (divergence of the inverse /// gradients): serves to calculate the phase space correction colvarvalue jd; }; inline colvarvalue::Type colvar::cvc::type() const { return x.type(); } inline colvarvalue colvar::cvc::value() const { return x; } inline colvarvalue colvar::cvc::system_force() const { return ft; } inline colvarvalue colvar::cvc::Jacobian_derivative() const { return jd; } inline cvm::real colvar::cvc::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { return x1.dist2 (x2); } inline colvarvalue colvar::cvc::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { return x1.dist2_grad (x2); } inline colvarvalue colvar::cvc::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { return x2.dist2_grad (x1); } inline cvm::real colvar::cvc::compare (colvarvalue const &x1, colvarvalue const &x2) const { if (this->type() == colvarvalue::type_scalar) { return cvm::real (x1 - x2); } else { cvm::fatal_error ("Error: you requested an operation which requires " "comparison between two non-scalar values.\n"); return 0.0; } } inline void colvar::cvc::wrap (colvarvalue &x) const { return; } /// \brief Colvar component: distance between the centers of mass of /// two groups (colvarvalue::type_scalar type, range [0:*)) /// /// This class also serves as the template for many collective /// variables with two atom groups: in this case, the /// distance::distance() constructor should be called on the same /// configuration string, to make the same member data and functions /// available to the derived object class colvar::distance : public colvar::cvc { protected: /// First atom group cvm::atom_group group1; /// Second atom group cvm::atom_group group2; /// Vector distance, cached to be recycled cvm::rvector dist_v; /// Use absolute positions, ignoring PBCs when present bool b_no_PBC; /// Compute system force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) bool b_1site_force; public: distance (std::string const &conf, bool twogroups = true); distance(); virtual inline ~distance() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; // \brief Colvar component: distance vector between centers of mass // of two groups (\link colvarvalue::type_vector \endlink type, // range (-*:*)x(-*:*)x(-*:*)) class colvar::distance_vec : public colvar::distance { public: distance_vec (std::string const &conf); distance_vec(); virtual inline ~distance_vec() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); /// Redefined to handle the box periodicity virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the box periodicity virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the box periodicity virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the box periodicity virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: distance unit vector (direction) between /// centers of mass of two groups (colvarvalue::type_unitvector type, /// range [-1:1]x[-1:1]x[-1:1]) class colvar::distance_dir : public colvar::distance { public: distance_dir (std::string const &conf); distance_dir(); virtual inline ~distance_dir() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: projection of the distance vector along /// an axis (colvarvalue::type_scalar type, range (-*:*)) class colvar::distance_z : public colvar::cvc { protected: /// Main atom group cvm::atom_group main; /// Reference atom group cvm::atom_group ref1; /// Optional, second ref atom group cvm::atom_group ref2; /// Use absolute positions, ignoring PBCs when present bool b_no_PBC; /// Compute system force on one site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) bool b_1site_force; /// Vector on which the distance vector is projected cvm::rvector axis; /// Norm of the axis cvm::real axis_norm; /// Vector distance, cached to be recycled cvm::rvector dist_v; /// Flag: using a fixed axis vector? bool fixed_axis; public: distance_z (std::string const &conf); distance_z(); virtual inline ~distance_z() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Redefined to make use of the user-provided period virtual void wrap (colvarvalue &x) const; }; /// \brief Colvar component: projection of the distance vector on a /// plane (colvarvalue::type_scalar type, range [0:*)) class colvar::distance_xy : public colvar::distance_z { protected: /// Components of the distance vector orthogonal to the axis cvm::rvector dist_v_ortho; /// Vector distances cvm::rvector v12, v13; public: distance_xy (std::string const &conf); distance_xy(); virtual inline ~distance_xy() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power, /// as in NMR refinements (colvarvalue::type_scalar type, range (0:*)) class colvar::distance_inv : public colvar::distance { protected: /// Components of the distance vector orthogonal to the axis int exponent; public: distance_inv (std::string const &conf); distance_inv(); virtual inline ~distance_inv() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: Radius of gyration of an atom group /// (colvarvalue::type_scalar type, range [0:*)) class colvar::gyration : public colvar::cvc { protected: /// Atoms involved cvm::atom_group atoms; public: /// Constructor gyration (std::string const &conf); gyration(); virtual inline ~gyration() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: moment of inertia of an atom group /// (colvarvalue::type_scalar type, range [0:*)) class colvar::inertia - : public colvar::cvc + : public colvar::gyration { -protected: - /// Atoms involved - cvm::atom_group atoms; public: /// Constructor inertia (std::string const &conf); inertia(); virtual inline ~inertia() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: moment of inertia of an atom group /// around a user-defined axis (colvarvalue::type_scalar type, range [0:*)) class colvar::inertia_z : public colvar::inertia { protected: /// Vector on which the inertia tensor is projected cvm::rvector axis; public: /// Constructor inertia_z (std::string const &conf); inertia_z(); virtual inline ~inertia_z() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: projection of 3N coordinates onto an /// eigenvector (colvarvalue::type_scalar type, range (-*:*)) class colvar::eigenvector : public colvar::cvc { protected: /// Atom group cvm::atom_group atoms; /// Reference coordinates std::vector ref_pos; - /// Eigenvector (of a normal or essential mode) + /// Geometric center of the reference coordinates + cvm::rvector ref_pos_center; + + /// Eigenvector (of a normal or essential mode): will always have zero center std::vector eigenvec; /// Inverse square norm of the eigenvector cvm::real eigenvec_invnorm2; public: /// Constructor eigenvector (std::string const &conf); virtual inline ~eigenvector() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: angle between the centers of mass of /// three groups (colvarvalue::type_scalar type, range [0:PI]) class colvar::angle : public colvar::cvc { protected: /// Atom group cvm::atom_group group1; /// Atom group cvm::atom_group group2; /// Atom group cvm::atom_group group3; /// Inter site vectors cvm::rvector r21, r23; /// Inter site vector norms cvm::real r21l, r23l; /// Derivatives wrt group centers of mass cvm::rvector dxdr1, dxdr3; /// Compute system force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) /// (or to allow dummy atoms) bool b_1site_force; public: /// Initialize by parsing the configuration angle (std::string const &conf); /// \brief Initialize the three groups after three atoms angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3); angle(); virtual inline ~angle() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: dihedral between the centers of mass of /// four groups (colvarvalue::type_scalar type, range [-PI:PI]) class colvar::dihedral : public colvar::cvc { protected: /// Atom group cvm::atom_group group1; /// Atom group cvm::atom_group group2; /// Atom group cvm::atom_group group3; /// Atom group cvm::atom_group group4; /// Inter site vectors cvm::rvector r12, r23, r34; /// \brief Compute system force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) bool b_1site_force; public: /// Initialize by parsing the configuration dihedral (std::string const &conf); /// \brief Initialize the four groups after four atoms dihedral (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4); dihedral(); virtual inline ~dihedral() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); /// Redefined to handle the 2*PI periodicity virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual void wrap (colvarvalue &x) const; }; /// \brief Colvar component: coordination number between two groups /// (colvarvalue::type_scalar type, range [0:N1*N2]) class colvar::coordnum : public colvar::distance { protected: /// \brief "Cutoff" for isotropic calculation (default) cvm::real r0; /// \brief "Cutoff vector" for anisotropic calculation cvm::rvector r0_vec; /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be /// used bool b_anisotropic; /// Integer exponent of the function numerator int en; /// Integer exponent of the function denominator int ed; /// \brief If true, group2 will be treated as a single atom /// (default: loop over all pairs of atoms in group1 and group2) bool b_group2_center_only; public: /// Constructor coordnum (std::string const &conf); coordnum(); virtual inline ~coordnum() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); template /// \brief Calculate a coordination number through the function /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the /// coordination number \param exp_num \i n exponent \param exp_den /// \i m exponent \param A1 atom \param A2 atom static cvm::real switching_function (cvm::real const &r0, int const &exp_num, int const &exp_den, cvm::atom &A1, cvm::atom &A2); template /// \brief Calculate a coordination number through the function /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec /// vector of different cutoffs in the three directions \param /// exp_num \i n exponent \param exp_den \i m exponent \param A1 /// atom \param A2 atom static cvm::real switching_function (cvm::rvector const &r0_vec, int const &exp_num, int const &exp_den, cvm::atom &A1, cvm::atom &A2); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: self-coordination number within a group /// (colvarvalue::type_scalar type, range [0:N*(N-1)/2]) class colvar::selfcoordnum : public colvar::distance { protected: /// \brief "Cutoff" for isotropic calculation (default) cvm::real r0; /// Integer exponent of the function numerator int en; /// Integer exponent of the function denominator int ed; public: /// Constructor selfcoordnum (std::string const &conf); selfcoordnum(); virtual inline ~selfcoordnum() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); template /// \brief Calculate a coordination number through the function /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the /// coordination number \param exp_num \i n exponent \param exp_den /// \i m exponent \param A1 atom \param A2 atom static cvm::real switching_function (cvm::real const &r0, int const &exp_num, int const &exp_den, cvm::atom &A1, cvm::atom &A2); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: hydrogen bond, defined as the product of /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol)) /// (colvarvalue::type_scalar type, range [0:1]) class colvar::h_bond : public colvar::cvc { protected: /// Atoms involved in the component cvm::atom acceptor, donor; /// \brief "Cutoff" distance between acceptor and donor cvm::real r0; /// Integer exponent of the function numerator int en; /// Integer exponent of the function denominator int ed; public: h_bond (std::string const &conf); /// Constructor for atoms already allocated h_bond (cvm::atom const &acceptor, cvm::atom const &donor, cvm::real r0, int en, int ed); h_bond(); virtual ~h_bond(); virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: alpha helix content of a contiguous /// segment of 5 or more residues, implemented as a sum of phi/psi /// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type, /// range [0:1]) // class colvar::alpha_dihedrals // : public colvar::cvc // { // protected: // /// Alpha-helical reference phi value // cvm::real phi_ref; // /// Alpha-helical reference psi value // cvm::real psi_ref; // /// List of phi dihedral angles // std::vector phi; // /// List of psi dihedral angles // std::vector psi; // /// List of hydrogen bonds // std::vector hb; // public: // alpha_dihedrals (std::string const &conf); // alpha_dihedrals(); // virtual inline ~alpha_dihedrals() {} // virtual void calc_value(); // virtual void calc_gradients(); // virtual void apply_force (colvarvalue const &force); // virtual cvm::real dist2 (colvarvalue const &x1, // colvarvalue const &x2) const; // virtual colvarvalue dist2_lgrad (colvarvalue const &x1, // colvarvalue const &x2) const; // virtual colvarvalue dist2_rgrad (colvarvalue const &x1, // colvarvalue const &x2) const; // virtual cvm::real compare (colvarvalue const &x1, // colvarvalue const &x2) const; // }; /// \brief Colvar component: alpha helix content of a contiguous /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca /// angles and hydrogen bonds (colvarvalue::type_scalar type, range /// [0:1]) class colvar::alpha_angles : public colvar::cvc { protected: /// Reference Calpha-Calpha angle (default: 88 degrees) cvm::real theta_ref; /// Tolerance on the Calpha-Calpha angle cvm::real theta_tol; /// List of Calpha-Calpha angles std::vector theta; /// List of hydrogen bonds std::vector hb; /// Contribution of the hb terms cvm::real hb_coeff; public: alpha_angles (std::string const &conf); alpha_angles(); virtual ~alpha_angles(); void calc_value(); void calc_gradients(); void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: dihedPC /// Projection of the config onto a dihedral principal component /// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007) /// Based on a set of 'dihedral' cvcs class colvar::dihedPC : public colvar::cvc { protected: std::vector theta; std::vector coeffs; public: dihedPC (std::string const &conf); dihedPC(); virtual ~dihedPC(); void calc_value(); void calc_gradients(); void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: orientation in space of an atom group, /// with respect to a set of reference coordinates /// (colvarvalue::type_quaternion type, range /// [-1:1]x[-1:1]x[-1:1]x[-1:1]) class colvar::orientation : public colvar::cvc { protected: /// Atom group cvm::atom_group atoms; /// Center of geometry of the group cvm::atom_pos atoms_cog; /// Reference coordinates std::vector ref_pos; /// Rotation object cvm::rotation rot; /// \brief This is used to remove jumps in the sign of the /// quaternion, which may be annoying in the colvars trajectory cvm::quaternion ref_quat; public: orientation (std::string const &conf); orientation(); virtual inline ~orientation() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: angle of rotation with respect to a set /// of reference coordinates (colvarvalue::type_scalar type, range /// [0:PI)) class colvar::orientation_angle : public colvar::orientation { public: orientation_angle (std::string const &conf); orientation_angle(); virtual inline ~orientation_angle() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: projection of the orientation vector onto /// a predefined axis (colvarvalue::type_scalar type, range [-1:1]) class colvar::tilt : public colvar::orientation { protected: cvm::rvector axis; public: tilt (std::string const &conf); tilt(); virtual inline ~tilt() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; /// \brief Colvar component: angle of rotation around a predefined /// axis (colvarvalue::type_scalar type, range [-PI:PI]) class colvar::spin_angle : public colvar::orientation { protected: cvm::rvector axis; public: spin_angle (std::string const &conf); spin_angle(); virtual inline ~spin_angle() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); /// Redefined to handle the 2*PI periodicity virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity virtual void wrap (colvarvalue &x) const; }; /// \brief Colvar component: root mean square deviation (RMSD) of a /// group with respect to a set of reference coordinates; uses \link /// colvar::orientation \endlink to calculate the rotation matrix /// (colvarvalue::type_scalar type, range [0:*)) class colvar::rmsd : public colvar::cvc { protected: /// Atom group cvm::atom_group atoms; /// Reference coordinates (for RMSD calculation only) std::vector ref_pos; public: /// Constructor rmsd (std::string const &conf); virtual inline ~rmsd() {} virtual void calc_value(); virtual void calc_gradients(); virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force (colvarvalue const &force); virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const; virtual cvm::real compare (colvarvalue const &x1, colvarvalue const &x2) const; }; - // metrics functions for cvc implementations with a periodicity inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return diff * diff; } inline colvarvalue colvar::dihedral::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return 2.0 * diff; } inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return (-2.0) * diff; } inline cvm::real colvar::dihedral::compare (colvarvalue const &x1, colvarvalue const &x2) const { return dist2_lgrad (x1, x2); } inline void colvar::dihedral::wrap (colvarvalue &x) const { if ((x.real_value - wrap_center) >= 180.0) { x.real_value -= 360.0; return; } if ((x.real_value - wrap_center) < -180.0) { x.real_value += 360.0; return; } return; } inline cvm::real colvar::spin_angle::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return diff * diff; } inline colvarvalue colvar::spin_angle::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return 2.0 * diff; } inline colvarvalue colvar::spin_angle::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); return (-2.0) * diff; } inline cvm::real colvar::spin_angle::compare (colvarvalue const &x1, colvarvalue const &x2) const { return dist2_lgrad (x1, x2); } inline void colvar::spin_angle::wrap (colvarvalue &x) const { if ((x.real_value - wrap_center) >= 180.0) { x.real_value -= 360.0; return; } if ((x.real_value - wrap_center) < -180.0) { x.real_value += 360.0; return; } return; } // simple definitions of the distance functions; these are useful only // for optimization (the type check performed in the default // colvarcomp functions is skipped) // definitions assuming the scalar type #define simple_scalar_dist_functions(TYPE) \ \ inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1, \ colvarvalue const &x2) const \ { \ return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \ } \ \ inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1, \ colvarvalue const &x2) const \ { \ return 2.0 * (x1.real_value - x2.real_value); \ } \ \ inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1, \ colvarvalue const &x2) const \ { \ return this->dist2_lgrad (x2, x1); \ } \ \ inline cvm::real colvar::TYPE::compare (colvarvalue const &x1, \ colvarvalue const &x2) const \ { \ return this->dist2_lgrad (x1, x2); \ } \ \ simple_scalar_dist_functions (distance) // NOTE: distance_z has explicit functions, see below simple_scalar_dist_functions (distance_xy) simple_scalar_dist_functions (distance_inv) simple_scalar_dist_functions (angle) simple_scalar_dist_functions (coordnum) simple_scalar_dist_functions (selfcoordnum) simple_scalar_dist_functions (h_bond) simple_scalar_dist_functions (gyration) simple_scalar_dist_functions (inertia) simple_scalar_dist_functions (inertia_z) simple_scalar_dist_functions (rmsd) simple_scalar_dist_functions (orientation_angle) simple_scalar_dist_functions (tilt) simple_scalar_dist_functions (eigenvector) // simple_scalar_dist_functions (alpha_dihedrals) simple_scalar_dist_functions (alpha_angles) simple_scalar_dist_functions (dihedPC) // Projected distance // Differences should always be wrapped around 0 (ignoring wrap_center) inline cvm::real colvar::distance_z::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; if (period != 0.0) { cvm::real shift = std::floor (diff/period + 0.5); diff -= shift * period; } return diff * diff; } inline colvarvalue colvar::distance_z::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; if (period != 0.0) { cvm::real shift = std::floor (diff/period + 0.5); diff -= shift * period; } return 2.0 * diff; } inline colvarvalue colvar::distance_z::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; if (period != 0.0) { cvm::real shift = std::floor (diff/period + 0.5); diff -= shift * period; } return (-2.0) * diff; } inline cvm::real colvar::distance_z::compare (colvarvalue const &x1, colvarvalue const &x2) const { return dist2_lgrad (x1, x2); } inline void colvar::distance_z::wrap (colvarvalue &x) const { if (! this->b_periodic) { // don't wrap if the period has not been set return; } cvm::real shift = std::floor ((x.real_value - wrap_center) / period + 0.5); x.real_value -= shift * period; return; } // distance between three dimensional vectors // // TODO apply PBC to distance_vec // Note: differences should be centered around (0, 0, 0)! inline cvm::real colvar::distance_vec::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { return cvm::position_dist2 (x1.rvector_value, x2.rvector_value); } inline colvarvalue colvar::distance_vec::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); } inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); } inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1, colvarvalue const &x2) const { cvm::fatal_error ("Error: cannot compare() two distance vectors.\n"); return 0.0; } inline cvm::real colvar::distance_dir::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { return (x1.rvector_value - x2.rvector_value).norm2(); } inline colvarvalue colvar::distance_dir::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { return colvarvalue ((x1.rvector_value - x2.rvector_value), colvarvalue::type_unitvector); } inline colvarvalue colvar::distance_dir::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { return colvarvalue ((x2.rvector_value - x1.rvector_value), colvarvalue::type_unitvector); } inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1, colvarvalue const &x2) const { cvm::fatal_error ("Error: cannot compare() two distance directions.\n"); return 0.0; } // distance between quaternions inline cvm::real colvar::orientation::dist2 (colvarvalue const &x1, colvarvalue const &x2) const { return x1.quaternion_value.dist2 (x2); } inline colvarvalue colvar::orientation::dist2_lgrad (colvarvalue const &x1, colvarvalue const &x2) const { return x1.quaternion_value.dist2_grad (x2); } inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1, colvarvalue const &x2) const { return x2.quaternion_value.dist2_grad (x1); } inline cvm::real colvar::orientation::compare (colvarvalue const &x1, colvarvalue const &x2) const { cvm::fatal_error ("Error: cannot compare() two quaternions.\n"); return 0.0; } #endif // Emacs // Local Variables: // mode: C++ // End: diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index bc77c5a4b..9519620b4 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -1,1198 +1,1284 @@ #include #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" #include "colvar.h" #include "colvarcomp.h" -/// \file cvc_distance.cpp \brief Collective variables -/// determining various type of distances between two groups - // "twogroup" flag defaults to true; set to false by selfCoordNum // (only distance-derived component based on only one group) colvar::distance::distance (std::string const &conf, bool twogroups) : cvc (conf) { function_type = "distance"; b_inverse_gradients = true; b_Jacobian_derivative = true; if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) { 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"); } 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); } colvar::distance::distance() : cvc () { function_type = "distance"; b_inverse_gradients = true; b_Jacobian_derivative = true; b_1site_force = false; x.type (colvarvalue::type_scalar); } void colvar::distance::calc_value() { - group1.reset_atoms_data(); - group2.reset_atoms_data(); - - group1.read_positions(); - group2.read_positions(); - if (b_no_PBC) { dist_v = group2.center_of_mass() - group1.center_of_mass(); } else { dist_v = cvm::position_distance (group1.center_of_mass(), group2.center_of_mass()); } x.real_value = dist_v.norm(); } void colvar::distance::calc_gradients() { cvm::rvector const u = dist_v.unit(); group1.set_weighted_gradient (-1.0 * u); group2.set_weighted_gradient ( u); } void colvar::distance::calc_force_invgrads() { group1.read_system_forces(); if ( b_1site_force ) { ft.real_value = -1.0 * (group1.system_force() * dist_v.unit()); } else { group2.read_system_forces(); ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit()); } } void colvar::distance::calc_Jacobian_derivative() { jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0; } void colvar::distance::apply_force (colvarvalue const &force) { if (!group1.noforce) group1.apply_colvar_force (force.real_value); if (!group2.noforce) group2.apply_colvar_force (force.real_value); } colvar::distance_vec::distance_vec (std::string const &conf) : distance (conf) { function_type = "distance_vec"; x.type (colvarvalue::type_vector); } colvar::distance_vec::distance_vec() : distance() { function_type = "distance_vec"; x.type (colvarvalue::type_vector); } void colvar::distance_vec::calc_value() { - group1.reset_atoms_data(); - group2.reset_atoms_data(); - - group1.read_positions(); - group2.read_positions(); - if (b_no_PBC) { x.rvector_value = group2.center_of_mass() - group1.center_of_mass(); } else { x.rvector_value = cvm::position_distance (group1.center_of_mass(), group2.center_of_mass()); } } void colvar::distance_vec::calc_gradients() { // gradients are not stored: a 3x3 matrix for each atom would be // needed to store just the identity matrix } void colvar::distance_vec::apply_force (colvarvalue const &force) { if (!group1.noforce) group1.apply_force (-1.0 * force.rvector_value); if (!group2.noforce) group2.apply_force ( force.rvector_value); } colvar::distance_z::distance_z (std::string const &conf) : cvc (conf) { function_type = "distance_z"; b_inverse_gradients = true; b_Jacobian_derivative = true; x.type (colvarvalue::type_scalar); // TODO detect PBC from MD engine (in simple cases) // and then update period in real time if (period != 0.0) b_periodic = true; if ((wrap_center != 0.0) && (period == 0.0)) { cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component," " but its period has not been set.\n"); } parse_group (conf, "main", main); parse_group (conf, "ref", ref1); atom_groups.push_back (&main); atom_groups.push_back (&ref1); // this group is optional parse_group (conf, "ref2", ref2, true); if (ref2.size()) { atom_groups.push_back (&ref2); cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); fixed_axis = false; if (key_lookup (conf, "axis")) cvm::log ("Warning: explicit axis definition will be ignored!"); } else { if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) cvm::fatal_error ("Axis vector is zero!"); if (axis.norm2() != 1.0) { axis = axis.unit(); cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n"); } } fixed_axis = true; } if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) { 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"); } } colvar::distance_z::distance_z() { function_type = "distance_z"; b_inverse_gradients = true; b_Jacobian_derivative = true; x.type (colvarvalue::type_scalar); } void colvar::distance_z::calc_value() { - main.reset_atoms_data(); - ref1.reset_atoms_data(); - - main.read_positions(); - ref1.read_positions(); - if (fixed_axis) { if (b_no_PBC) { dist_v = main.center_of_mass() - ref1.center_of_mass(); } else { dist_v = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass()); } } else { - ref2.reset_atoms_data(); - ref2.read_positions(); if (b_no_PBC) { dist_v = main.center_of_mass() - (0.5 * (ref1.center_of_mass() + ref2.center_of_mass())); axis = ref2.center_of_mass() - ref1.center_of_mass(); } else { dist_v = cvm::position_distance (0.5 * (ref1.center_of_mass() + ref2.center_of_mass()), main.center_of_mass()); axis = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass()); } axis_norm = axis.norm(); axis = axis.unit(); } x.real_value = axis * dist_v; this->wrap (x); } void colvar::distance_z::calc_gradients() { main.set_weighted_gradient ( axis ); if (fixed_axis) { ref1.set_weighted_gradient (-1.0 * axis); } else { if (b_no_PBC) { ref1.set_weighted_gradient ( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() - x.real_value * axis )); ref2.set_weighted_gradient ( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() + x.real_value * axis )); } else { ref1.set_weighted_gradient ( 1.0 / axis_norm * ( cvm::position_distance (ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis )); ref2.set_weighted_gradient ( 1.0 / axis_norm * ( cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis )); } } } void colvar::distance_z::calc_force_invgrads() { main.read_system_forces(); if (fixed_axis && !b_1site_force) { ref1.read_system_forces(); ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis); } else { ft.real_value = main.system_force() * axis; } } void colvar::distance_z::calc_Jacobian_derivative() { jd.real_value = 0.0; } void colvar::distance_z::apply_force (colvarvalue const &force) { if (!ref1.noforce) ref1.apply_colvar_force (force.real_value); if (ref2.size() && !ref2.noforce) ref2.apply_colvar_force (force.real_value); if (!main.noforce) main.apply_colvar_force (force.real_value); } colvar::distance_xy::distance_xy (std::string const &conf) : distance_z (conf) { function_type = "distance_xy"; b_inverse_gradients = true; b_Jacobian_derivative = true; x.type (colvarvalue::type_scalar); } colvar::distance_xy::distance_xy() : distance_z() { function_type = "distance_xy"; b_inverse_gradients = true; b_Jacobian_derivative = true; x.type (colvarvalue::type_scalar); } void colvar::distance_xy::calc_value() { - ref1.reset_atoms_data(); - main.reset_atoms_data(); - - ref1.read_positions(); - main.read_positions(); - if (b_no_PBC) { dist_v = main.center_of_mass() - ref1.center_of_mass(); } else { dist_v = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass()); } if (!fixed_axis) { - ref2.reset_atoms_data(); - ref2.read_positions(); - if (b_no_PBC) { v12 = ref2.center_of_mass() - ref1.center_of_mass(); } else { v12 = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass()); } axis_norm = v12.norm(); axis = v12.unit(); } dist_v_ortho = dist_v - (dist_v * axis) * axis; x.real_value = dist_v_ortho.norm(); } void colvar::distance_xy::calc_gradients() { // Intermediate quantity (r_P3 / r_12 where P is the projection // of 3 (main) on the plane orthogonal to 12, containing 1 (ref1)) cvm::real A; cvm::real x_inv; if (x.real_value == 0.0) return; x_inv = 1.0 / x.real_value; if (fixed_axis) { ref1.set_weighted_gradient (-1.0 * x_inv * dist_v_ortho); main.set_weighted_gradient ( x_inv * dist_v_ortho); } else { if (b_no_PBC) { v13 = main.center_of_mass() - ref1.center_of_mass(); } else { v13 = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass()); } A = (dist_v * axis) / axis_norm; ref1.set_weighted_gradient ( (A - 1.0) * 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); } } void colvar::distance_xy::calc_force_invgrads() { main.read_system_forces(); if (fixed_axis && !b_1site_force) { ref1.read_system_forces(); ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho); } else { ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho; } } void colvar::distance_xy::calc_Jacobian_derivative() { jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0; } void colvar::distance_xy::apply_force (colvarvalue const &force) { if (!ref1.noforce) ref1.apply_colvar_force (force.real_value); if (ref2.size() && !ref2.noforce) ref2.apply_colvar_force (force.real_value); if (!main.noforce) main.apply_colvar_force (force.real_value); } colvar::distance_dir::distance_dir (std::string const &conf) : distance (conf) { function_type = "distance_dir"; x.type (colvarvalue::type_unitvector); } colvar::distance_dir::distance_dir() : distance() { function_type = "distance_dir"; x.type (colvarvalue::type_unitvector); } void colvar::distance_dir::calc_value() { - group1.reset_atoms_data(); - group2.reset_atoms_data(); - - group1.read_positions(); - group2.read_positions(); - if (b_no_PBC) { dist_v = group2.center_of_mass() - group1.center_of_mass(); } else { dist_v = cvm::position_distance (group1.center_of_mass(), group2.center_of_mass()); } x.rvector_value = dist_v.unit(); } void colvar::distance_dir::calc_gradients() { // gradients are computed on the fly within apply_force() // Note: could be a problem if a future bias relies on gradient // calculations... } void colvar::distance_dir::apply_force (colvarvalue const &force) { // remove the radial force component cvm::real const iprod = force.rvector_value * x.rvector_value; cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value; if (!group1.noforce) group1.apply_force (-1.0 * force_tang); if (!group2.noforce) group2.apply_force ( force_tang); } colvar::distance_inv::distance_inv (std::string const &conf) : distance (conf) { function_type = "distance_inv"; get_keyval (conf, "exponent", exponent, 6); if (exponent%2) { cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n"); } if (exponent <= 0) { cvm::fatal_error ("Error: negative or zero exponent provided.\n"); } for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { if (ai1->id == ai2->id) cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n"); } } b_inverse_gradients = false; b_Jacobian_derivative = false; x.type (colvarvalue::type_scalar); } colvar::distance_inv::distance_inv() { function_type = "distance_inv"; exponent = 6; b_inverse_gradients = false; b_Jacobian_derivative = false; b_1site_force = false; x.type (colvarvalue::type_scalar); } void colvar::distance_inv::calc_value() { - group1.reset_atoms_data(); - group2.reset_atoms_data(); - - group1.read_positions(); - group2.read_positions(); - x.real_value = 0.0; if (b_no_PBC) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { cvm::rvector const dv = ai2->pos - ai1->pos; cvm::real const d2 = dv.norm2(); cvm::real dinv = 1.0; for (int ne = 0; ne < exponent/2; ne++) dinv *= 1.0/d2; x.real_value += dinv; cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv; ai1->grad += -1.0 * dsumddv; ai2->grad += dsumddv; } } } else { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos); cvm::real const d2 = dv.norm2(); cvm::real dinv = 1.0; for (int ne = 0; ne < exponent/2; ne++) dinv *= 1.0/d2; x.real_value += dinv; cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv; ai1->grad += -1.0 * dsumddv; ai2->grad += dsumddv; } } } x.real_value *= 1.0 / cvm::real (group1.size() * group2.size()); x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent))); } 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()); for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { ai1->grad *= dxdsum; } for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { ai2->grad *= dxdsum; } } void colvar::distance_inv::apply_force (colvarvalue const &force) { if (!group1.noforce) group1.apply_colvar_force (force.real_value); if (!group2.noforce) group2.apply_colvar_force (force.real_value); } colvar::gyration::gyration (std::string const &conf) : cvc (conf) { function_type = "gyration"; b_inverse_gradients = true; b_Jacobian_derivative = true; parse_group (conf, "atoms", atoms); atom_groups.push_back (&atoms); + + if (atoms.b_user_defined_fit) { + cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); + } else { + // default: fit everything + cvm::log ("No explicit fitting parameters: enabling centerReference by default to align with the origin\n"); + cvm::log (" (rotateReference will not affect results for gyration and inertia, but it will for inertiaZ).\n"); + atoms.b_center = true; + atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0)); + } + x.type (colvarvalue::type_scalar); } colvar::gyration::gyration() { function_type = "gyration"; b_inverse_gradients = true; b_Jacobian_derivative = true; x.type (colvarvalue::type_scalar); } void colvar::gyration::calc_value() { - atoms.reset_atoms_data(); - atoms.read_positions(); - atoms.apply_translation ((-1.0) * atoms.center_of_geometry()); - x.real_value = 0.0; for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { x.real_value += (ai->pos).norm2(); } x.real_value = std::sqrt (x.real_value / cvm::real (atoms.size())); } void colvar::gyration::calc_gradients() { cvm::real const drdx = 1.0/(cvm::real (atoms.size()) * x.real_value); for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { ai->grad = drdx * ai->pos; } + + if (b_debug_gradients) { + cvm::log ("Debugging gradients:\n"); + debug_gradients (atoms); + } } void colvar::gyration::calc_force_invgrads() { atoms.read_system_forces(); cvm::real const dxdr = 1.0/x.real_value; ft.real_value = 0.0; for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { ft.real_value += dxdr * ai->pos * ai->system_force; } } void colvar::gyration::calc_Jacobian_derivative() { 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) { if (!atoms.noforce) atoms.apply_colvar_force (force.real_value); } colvar::inertia::inertia (std::string const &conf) - : cvc (conf) + : gyration (conf) { function_type = "inertia"; - parse_group (conf, "atoms", atoms); - atom_groups.push_back (&atoms); + b_inverse_gradients = false; + b_Jacobian_derivative = false; x.type (colvarvalue::type_scalar); } colvar::inertia::inertia() { function_type = "inertia"; x.type (colvarvalue::type_scalar); } void colvar::inertia::calc_value() { - atoms.reset_atoms_data(); - atoms.read_positions(); - atoms.apply_translation ((-1.0) * atoms.center_of_mass()); - x.real_value = 0.0; for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { x.real_value += ai->mass * (ai->pos).norm2(); } } void colvar::inertia::calc_gradients() { for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { ai->grad = 2.0 * ai->mass * ai->pos; } + + if (b_debug_gradients) { + cvm::log ("Debugging gradients:\n"); + debug_gradients (atoms); + } } void colvar::inertia::apply_force (colvarvalue const &force) { if (!atoms.noforce) atoms.apply_colvar_force (force.real_value); } colvar::inertia_z::inertia_z (std::string const &conf) : inertia (conf) { function_type = "inertia_z"; if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) cvm::fatal_error ("Axis vector is zero!"); if (axis.norm2() != 1.0) { axis = axis.unit(); cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n"); } } x.type (colvarvalue::type_scalar); } colvar::inertia_z::inertia_z() { function_type = "inertia_z"; x.type (colvarvalue::type_scalar); } void colvar::inertia_z::calc_value() { - atoms.reset_atoms_data(); - atoms.read_positions(); - atoms.apply_translation ((-1.0) * atoms.center_of_mass()); - x.real_value = 0.0; for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { cvm::real const iprod = ai->pos * axis; x.real_value += ai->mass * iprod * iprod; } } void colvar::inertia_z::calc_gradients() { for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { ai->grad = 2.0 * ai->mass * (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) { if (!atoms.noforce) atoms.apply_colvar_force (force.real_value); } colvar::rmsd::rmsd (std::string const &conf) : cvc (conf) { b_inverse_gradients = true; b_Jacobian_derivative = true; function_type = "rmsd"; x.type (colvarvalue::type_scalar); parse_group (conf, "atoms", atoms); atom_groups.push_back (&atoms); if (atoms.b_dummy) cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom."); - if (atoms.b_user_defined_fit) { - cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); - cvm::log ("Not computing standard minimum RMSD (if that is the desired result, leave atom group parameters at default values)."); - } else { - // Default: fit everything - // NOTE: this won't work for class V cvc's - atoms.b_center = true; - atoms.b_rotate = true; - } - - if (atoms.b_center || atoms.b_rotate) { - // request the calculation of the derivatives of the rotation defined by the atom group - atoms.rot.request_group1_gradients (atoms.size()); - // request derivatives of optimal rotation wrt 2nd group for Jacobian: this is only - // required for ABF, but we do both groups here for better caching - atoms.rot.request_group2_gradients (atoms.size()); - - if (atoms.ref_pos_group != NULL) { - cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " - "Jacobian derivatives of the RMSD will not be calculated.\n"); - b_Jacobian_derivative = false; - } + if (atoms.ref_pos_group != NULL) { + cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " + "Jacobian derivatives of the RMSD will not be calculated.\n"); + b_Jacobian_derivative = false; } - // the following is a simplified version of the corresponding atom group options: - // to define the reference coordinates to compute this variable + // the following is a simplified version of the corresponding atom group options; + // 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 if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) { cvm::log ("Using reference positions from configuration file to calculate the variable.\n"); if (ref_pos.size() != atoms.size()) { cvm::fatal_error ("Error: the number of reference positions provided ("+ cvm::to_str (ref_pos.size())+ ") does not match the number of atoms of group \"atoms\" ("+ cvm::to_str (atoms.size())+").\n"); } } { std::string ref_pos_file; if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) { if (ref_pos.size()) { cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and " "\"refPositions\" at the same time.\n"); } std::string ref_pos_col; double ref_pos_col_value; if (get_keyval (conf, "refPositionsCol", ref_pos_col, std::string (""))) { // if provided, use PDB column to select coordinates bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0); if (found && !ref_pos_col_value) cvm::fatal_error ("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); } else { // if not, rely on existing atom indices for the group atoms.create_sorted_ids(); } ref_pos.resize (atoms.size()); cvm::load_coords (ref_pos_file.c_str(), ref_pos, atoms.sorted_ids, ref_pos_col, ref_pos_col_value); } } - if (!atoms.b_user_defined_fit) { - // Default case: reference positions for calculating the rmsd are also those used + if (atoms.b_user_defined_fit) { + cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); + } else { + // Default: fit everything + cvm::log ("No explicit fitting parameters: enabling both centerReference and rotateReference to compute standard minimum RMSD."); + atoms.b_center = true; + atoms.b_rotate = true; + // default case: reference positions for calculating the rmsd are also those used // for fitting atoms.ref_pos = ref_pos; atoms.center_ref_pos(); + + // request the calculation of the derivatives of the rotation defined by the atom group + atoms.rot.request_group1_gradients (atoms.size()); + // 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 + atoms.rot.request_group2_gradients (atoms.size()); } } void colvar::rmsd::calc_value() { - atoms.reset_atoms_data(); - atoms.read_positions(); - // rotational fit is done internally - - // atoms_cog = atoms.center_of_geometry(); - // rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog)); - - // cvm::real group_pos_sum2 = 0.0; - // for (size_t i = 0; i < atoms.size(); i++) { - // group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2(); - // } - - // // value of the RMSD (Coutsias et al) - // cvm::real const MSD = 1.0/(cvm::real (atoms.size())) * - // ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda ); - - // x.real_value = (MSD > 0.0) ? std::sqrt (MSD) : 0.0; + // rotational-translational fit is handled by the atom group x.real_value = 0.0; for (size_t ia = 0; ia < atoms.size(); ia++) { x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2(); } x.real_value /= cvm::real (atoms.size()); // MSD - x.real_value = (x.real_value > 0.0) ? std::sqrt (x.real_value) : 0.0; + x.real_value = std::sqrt (x.real_value); } void colvar::rmsd::calc_gradients() { cvm::real const drmsddx2 = (x.real_value > 0.0) ? 0.5 / (x.real_value * cvm::real (atoms.size())) : 0.0; for (size_t ia = 0; ia < atoms.size(); 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) { if (!atoms.noforce) atoms.apply_colvar_force (force.real_value); } void colvar::rmsd::calc_force_invgrads() { atoms.read_system_forces(); ft.real_value = 0.0; // Note: gradient square norm is 1/N_atoms for (size_t ia = 0; ia < atoms.size(); ia++) { ft.real_value += atoms[ia].grad * atoms[ia].system_force; } ft.real_value *= atoms.size(); } void colvar::rmsd::calc_Jacobian_derivative() { // divergence of the rotated coordinates (including only derivatives of the rotation matrix) cvm::real divergence = 0.0; if (atoms.b_rotate) { // gradient of the rotation matrix cvm::matrix2d grad_rot_mat; // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; for (size_t ia = 0; ia < atoms.size(); ia++) { // Gradient of optimal quaternion wrt current Cartesian position cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia]; g11 = 2.0 * (atoms.rot.q)[1]*dq[1]; g22 = 2.0 * (atoms.rot.q)[2]*dq[2]; g33 = 2.0 * (atoms.rot.q)[3]*dq[3]; 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]; 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]; 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]; // Gradient of the rotation matrix wrt current Cartesian position grad_rot_mat[0][0] = -2.0 * (g22 + g33); grad_rot_mat[1][0] = 2.0 * (g12 + g03); grad_rot_mat[2][0] = 2.0 * (g13 - g02); grad_rot_mat[0][1] = 2.0 * (g12 - g03); grad_rot_mat[1][1] = -2.0 * (g11 + g33); grad_rot_mat[2][1] = 2.0 * (g01 + g23); grad_rot_mat[0][2] = 2.0 * (g02 + g13); grad_rot_mat[1][2] = 2.0 * (g23 - g01); grad_rot_mat[2][2] = -2.0 * (g11 + g22); - cvm::atom_pos &x = atoms[ia].pos; + cvm::atom_pos &y = ref_pos[ia]; - for (size_t i = 0; i < 3; i++) { - for (size_t j = 0; j < 3; j++) { - divergence += grad_rot_mat[i][j][i] * x[j]; + for (size_t alpha = 0; alpha < 3; alpha++) { + for (size_t beta = 0; beta < 3; beta++) { + divergence += grad_rot_mat[beta][alpha][alpha] * y[beta]; + // Note: equation was derived for inverse rotation (see colvars paper) + // so here the matrix is transposed + // (eq would give divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];) } } } } jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0; } colvar::eigenvector::eigenvector (std::string const &conf) : cvc (conf) { b_inverse_gradients = true; b_Jacobian_derivative = true; function_type = "eigenvector"; x.type (colvarvalue::type_scalar); parse_group (conf, "atoms", atoms); atom_groups.push_back (&atoms); - if (atoms.b_rotate) { - cvm::fatal_error ("Error: rotateReference should be disabled:" - "eigenvector component will set it internally."); - } - if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) { - cvm::log ("Using reference positions from input file.\n"); - if (ref_pos.size() != atoms.size()) { - cvm::fatal_error ("Error: reference positions do not " - "match the number of requested atoms.\n"); + { + bool const b_inline = get_keyval (conf, "refPositions", ref_pos, ref_pos); + + if (b_inline) { + cvm::log ("Using reference positions from input file.\n"); + if (ref_pos.size() != atoms.size()) { + cvm::fatal_error ("Error: reference positions do not " + "match the number of requested atoms.\n"); + } } - } - { std::string file_name; if (get_keyval (conf, "refPositionsFile", file_name)) { + if (b_inline) + cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n"); + std::string file_col; double file_col_value; if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) { // use PDB flags if column is provided bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0); if (found && !file_col_value) cvm::fatal_error ("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); } else { // if not, use atom indices atoms.create_sorted_ids(); } ref_pos.resize (atoms.size()); cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value); } } - // Set mobile frame of reference for atom group - atoms.b_center = true; - atoms.b_rotate = true; - atoms.ref_pos = ref_pos; - atoms.center_ref_pos(); - - // now load the eigenvector - if (get_keyval (conf, "vector", eigenvec, eigenvec)) { - cvm::log ("Using vector components from input file.\n"); - if (eigenvec.size() != atoms.size()) { - cvm::fatal_error ("Error: vector components do not " - "match the number of requested atoms.\n"); - } + // 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); + for (size_t i = 0; i < atoms.size(); i++) { + ref_pos_center += ref_pos[i]; + } + ref_pos_center *= 1.0 / atoms.size(); + + if (atoms.b_user_defined_fit) { + cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); + } else { + // default: fit everything + cvm::log ("Enabling centerReference and rotateReference (minimize RMSD before calculating the projection).\n"); + atoms.b_center = true; + atoms.b_rotate = true; + atoms.ref_pos = ref_pos; + atoms.center_ref_pos(); + + // request the calculation of the derivatives of the rotation defined by the atom group + atoms.rot.request_group1_gradients (atoms.size()); + // 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 + atoms.rot.request_group2_gradients (atoms.size()); } { + bool const b_inline = get_keyval (conf, "vector", eigenvec, eigenvec); + // now load the eigenvector + if (b_inline) { + cvm::log ("Using vector components from input file.\n"); + if (eigenvec.size() != atoms.size()) { + cvm::fatal_error ("Error: vector components do not " + "match the number of requested atoms.\n"); + } + } + std::string file_name; if (get_keyval (conf, "vectorFile", file_name)) { - std::string file_col; - if (!get_keyval (conf, "vectorCol", file_col, std::string (""))) { - cvm::fatal_error ("Error: parameter vectorCol is required if vectorFile is set.\n"); - } + if (b_inline) + cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n"); + std::string file_col; double file_col_value; - bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0); - if (found && !file_col_value) - cvm::fatal_error ("Error: eigenvectorColValue, " - "if provided, must be non-zero.\n"); + if (get_keyval (conf, "vectorCol", file_col, std::string (""))) { + // use PDB flags if column is provided + bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0); + if (found && !file_col_value) + cvm::fatal_error ("Error: vectorColValue, " + "if provided, must be non-zero.\n"); + } else { + // if not, use atom indices + atoms.create_sorted_ids(); + } eigenvec.resize (atoms.size()); cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value); } } if (!ref_pos.size() || !eigenvec.size()) { cvm::fatal_error ("Error: both reference coordinates" "and eigenvector must be defined.\n"); } - cvm::rvector center (0.0, 0.0, 0.0); - eigenvec_invnorm2 = 0.0; - + cvm::rvector eig_center (0.0, 0.0, 0.0); for (size_t i = 0; i < atoms.size(); i++) { - center += eigenvec[i]; + eig_center += eigenvec[i]; + } + eig_center *= 1.0 / atoms.size(); + cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n"); + + bool b_vector_difference; + get_keyval (conf, "vectorDifference", b_vector_difference, false); + if (b_vector_difference) { + + cvm::log ("Subtracting the reference positions from the provided vector.\n"); + + if (atoms.b_center) { + cvm::log ("centerReference is on: recentering the provided vector to the reference positions.\n"); + // both sets should be centered on the origin for fitting + for (size_t i = 0; i < atoms.size(); i++) { + eigenvec[i] -= eig_center; + ref_pos[i] -= ref_pos_center; + } + } + + if (atoms.b_rotate) { + cvm::log ("rotateReference is on: aligning the provided vector to the reference positions.\n"); + atoms.rot.calc_optimal_rotation (eigenvec, ref_pos); + for (size_t i = 0; i < atoms.size(); i++) { + eigenvec[i] = atoms.rot.rotate (eigenvec[i]); + eigenvec[i] -= ref_pos[i]; + } + } + + if (atoms.b_center) { + // bring back the ref positions to where they were + for (size_t i = 0; i < atoms.size(); i++) { + ref_pos[i] -= ref_pos_center; + } + } + + } else { + cvm::log ("Centering the provided vector to zero.\n"); + for (size_t i = 0; i < atoms.size(); i++) { + eigenvec[i] -= eig_center; + } } - center *= 1.0 / atoms.size(); - cvm::log ("Subtracting center of eigenvector components: " + cvm::to_str (center) + "\n"); + cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n"); + // for inverse gradients + eigenvec_invnorm2 = 0.0; for (size_t i = 0; i < atoms.size(); i++) { - eigenvec[i] = eigenvec[i] - center; eigenvec_invnorm2 += eigenvec[i].norm2(); } eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2; - - // request derivatives of optimal rotation wrt 2nd group - // for Jacobian - atoms.rot.request_group1_gradients(atoms.size()); - atoms.rot.request_group2_gradients(atoms.size()); } void colvar::eigenvector::calc_value() { - atoms.reset_atoms_data(); - atoms.read_positions(); // this will also update atoms.rot - x.real_value = 0.0; for (size_t i = 0; i < atoms.size(); i++) { x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i]; } } void colvar::eigenvector::calc_gradients() { // There are two versions of this code // The simple version is not formally exact, but its // results are numerically indistinguishable from the // exact one. The exact one is more expensive and possibly // less stable in cases where the optimal rotation // becomes ill-defined. // Version A: simple, intuitive, cheap, robust. Wrong. // Works just fine in practice. for (size_t ia = 0; ia < atoms.size(); ia++) { atoms[ia].grad = eigenvec[ia]; } -/* + std::vector gradients (atoms.size()); + for (size_t i = 0; i < atoms.size(); i++) { + gradients[i] = atoms[i].grad; + } + std::vector gradients_jh_correct (atoms.size()); + + /* + // Version B: complex, expensive, fragile. Right. // gradient of the rotation matrix cvm::matrix2d grad_rot_mat; cvm::quaternion &quat0 = atoms.rot.q; // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; // a term that involves the rotation gradients cvm::rvector rot_grad_term; cvm::atom_pos x_relative; cvm::atom_pos atoms_cog = atoms.center_of_geometry(); for (size_t ia = 0; ia < atoms.size(); ia++) { // Gradient of optimal quaternion wrt current Cartesian position - // WARNING: we want derivatives wrt the FIRST group here (unlike RMSD) cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia]; g11 = 2.0 * quat0[1]*dq[1]; g22 = 2.0 * quat0[2]*dq[2]; g33 = 2.0 * quat0[3]*dq[3]; g01 = quat0[0]*dq[1] + quat0[1]*dq[0]; g02 = quat0[0]*dq[2] + quat0[2]*dq[0]; g03 = quat0[0]*dq[3] + quat0[3]*dq[0]; g12 = quat0[1]*dq[2] + quat0[2]*dq[1]; g13 = quat0[1]*dq[3] + quat0[3]*dq[1]; g23 = quat0[2]*dq[3] + quat0[3]*dq[2]; // Gradient of the rotation matrix wrt current Cartesian position grad_rot_mat[0][0] = -2.0 * (g22 + g33); grad_rot_mat[1][0] = 2.0 * (g12 + g03); grad_rot_mat[2][0] = 2.0 * (g13 - g02); grad_rot_mat[0][1] = 2.0 * (g12 - g03); grad_rot_mat[1][1] = -2.0 * (g11 + g33); grad_rot_mat[2][1] = 2.0 * (g01 + g23); grad_rot_mat[0][2] = 2.0 * (g02 + g13); grad_rot_mat[1][2] = 2.0 * (g23 - g01); grad_rot_mat[2][2] = -2.0 * (g11 + g22); // this term needs to be rotated back, so we sum it separately rot_grad_term.reset(); for (size_t ja = 0; ja < atoms.size(); ja++) { x_relative = atoms[ja].pos - atoms_cog; for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { rot_grad_term += eigenvec[ja][i] * grad_rot_mat[i][j] * x_relative[j]; } } } // Rotate correction term back to reference frame - atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term); + // atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term); + // atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term); + + gradients_jh_correct[ia] = eigenvec[ia] + quat0.rotate (rot_grad_term); + } + */ + + // TODO replace this with a call to debug_gradients() + if (b_debug_gradients) { + + // the following code should work for any scalar variable; + // the only difference is the name of the atom group object (here, "atoms") + + cvm::rotation const rot_0 = atoms.rot; + cvm::rotation const rot_inv = atoms.rot.inverse(); + + cvm::real const x_0 = x.real_value; + + cvm::log ("gradients = "+cvm::to_str (gradients)+"\n"); + /* cvm::log ("gradients(JH-correct) = "+cvm::to_str (gradients_jh_correct)+"\n"); */ + + if (atoms.b_fit_gradients) { + atoms.calc_fit_gradients(); + if (atoms.b_rotate) { + // fit_gradients are in the original frame, output them in the rotated frame + for (size_t j = 0; j < atoms.fit_gradients.size(); j++) { + atoms.fit_gradients[j] = rot_0.rotate (atoms.fit_gradients[j]); + } + } + cvm::log ("fit_gradients = "+cvm::to_str (atoms.fit_gradients)+"\n"); + if (atoms.b_rotate) { + for (size_t j = 0; j < atoms.fit_gradients.size(); j++) { + atoms.fit_gradients[j] = rot_inv.rotate (atoms.fit_gradients[j]); + } + } + } + + for (size_t ia = 0; ia < atoms.size(); ia++) { + + // tests are best conducted in the unrotated (simulation) frame + cvm::rvector const atom_grad = atoms.b_rotate ? + rot_inv.rotate (atoms[ia].grad) : + atoms[ia].grad; + + for (size_t id = 0; id < 3; id++) { + // (re)read original positions + atoms.read_positions(); + // change one coordinate + atoms[ia].pos[id] += cvm::debug_gradients_step_size; + // (re)do the fit (if defined) + if (atoms.b_center || atoms.b_rotate) { + atoms.calc_apply_roto_translation(); + } + calc_value(); + cvm::real const x_1 = x.real_value; + cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n"); + cvm::log ("dx(real) = "+cvm::to_str (x_1 - x_0)+"\n"); + cvm::real const dx_pred = atoms.b_fit_gradients ? + (cvm::debug_gradients_step_size * (atom_grad[id] + atoms.fit_gradients[ia][id])) : + (cvm::debug_gradients_step_size * atom_grad[id]); + cvm::log ("dx(pred) = "+cvm::to_str (dx_pred)+"\n"); + cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+ + cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) / + std::fabs (x_1 - x_0), + cvm::cv_width, cvm::cv_prec)+ + ".\n"); + + } + } } -*/ } void colvar::eigenvector::apply_force (colvarvalue const &force) { if (!atoms.noforce) atoms.apply_colvar_force (force.real_value); } void colvar::eigenvector::calc_force_invgrads() { atoms.read_system_forces(); ft.real_value = 0.0; for (size_t ia = 0; ia < atoms.size(); ia++) { ft.real_value += eigenvec_invnorm2 * atoms[ia].grad * atoms[ia].system_force; } } void colvar::eigenvector::calc_Jacobian_derivative() { // gradient of the rotation matrix cvm::matrix2d grad_rot_mat; cvm::quaternion &quat0 = atoms.rot.q; // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; cvm::atom_pos x_relative; cvm::real sum = 0.0; for (size_t ia = 0; ia < atoms.size(); ia++) { // Gradient of optimal quaternion wrt current Cartesian position // trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t // we can just transpose the derivatives of the direct matrix cvm::vector1d< cvm::rvector, 4 > &dq_1 = atoms.rot.dQ0_1[ia]; g11 = 2.0 * quat0[1]*dq_1[1]; g22 = 2.0 * quat0[2]*dq_1[2]; g33 = 2.0 * quat0[3]*dq_1[3]; g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0]; g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0]; g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0]; g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1]; g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1]; g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2]; // Gradient of the inverse rotation matrix wrt current Cartesian position // (transpose of the gradient of the direct rotation) grad_rot_mat[0][0] = -2.0 * (g22 + g33); grad_rot_mat[0][1] = 2.0 * (g12 + g03); grad_rot_mat[0][2] = 2.0 * (g13 - g02); grad_rot_mat[1][0] = 2.0 * (g12 - g03); grad_rot_mat[1][1] = -2.0 * (g11 + g33); grad_rot_mat[1][2] = 2.0 * (g01 + g23); grad_rot_mat[2][0] = 2.0 * (g02 + g13); grad_rot_mat[2][1] = 2.0 * (g23 - g01); grad_rot_mat[2][2] = -2.0 * (g11 + g22); for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { sum += grad_rot_mat[i][j][i] * eigenvec[ia][j]; } } } jd.real_value = sum * std::sqrt (eigenvec_invnorm2); } + diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index c6e4d780a..80a7b70b1 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -1,1187 +1,1188 @@ // -*- c++ -*- #ifndef COLVARGRID_H #define COLVARGRID_H #include #include #include #include "colvar.h" #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" /// \brief Grid of values of a function of several collective /// variables \param T The data type /// /// Only scalar colvars supported so far template class colvar_grid : public colvarparse { protected: /// Number of dimensions size_t nd; /// Number of points along each dimension std::vector nx; /// Cumulative number of points along each dimension std::vector nxc; /// \brief Multiplicity of each datum (allow the binning of /// non-scalar types) size_t mult; /// Total number of grid points size_t nt; /// Low-level array of values std::vector data; /// Newly read data (used for count grids, when adding several grids read from disk) std::vector new_data; /// Colvars collected in this grid std::vector cv; /// Do we request actual value (for extended-system colvars)? std::vector actual_value; /// Get the low-level index corresponding to an index inline size_t address (std::vector const &ix) const { size_t addr = 0; for (size_t i = 0; i < nd; i++) { addr += ix[i]*nxc[i]; if (cvm::debug()) { if (ix[i] >= nx[i]) cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n"); } } return addr; } public: /// Lower boundaries of the colvars in this grid std::vector lower_boundaries; /// Upper boundaries of the colvars in this grid std::vector upper_boundaries; /// Whether some colvars are periodic std::vector periodic; /// Whether some colvars have hard lower boundaries std::vector hard_lower_boundaries; /// Whether some colvars have hard upper boundaries std::vector hard_upper_boundaries; /// Widths of the colvars in this grid std::vector widths; /// True if this is a count grid related to another grid of data bool has_parent_data; /// Whether this grid has been filled with data or is still empty bool has_data; /// Return the number of colvars inline size_t number_of_colvars() const { return nd; } /// Return the number of points in the i-th direction, if provided, or /// the total number inline size_t number_of_points (int const icv = -1) const { if (icv < 0) { return nt; } else { return nx[icv]; } } /// Get the sizes in each direction inline std::vector const & sizes() const { return nx; } /// Set the sizes in each direction inline void set_sizes (std::vector const &new_sizes) { nx = new_sizes; } /// Return the multiplicity of the type used inline size_t multiplicity() const { return mult; } /// \brief Allocate data (allow initialization also after construction) void create (std::vector const &nx_i, T const &t = T(), size_t const &mult_i = 1) { mult = mult_i; nd = nx_i.size(); nxc.resize (nd); nx = nx_i; nt = mult; for (int i = nd-1; i >= 0; i--) { if (nx_i[i] <= 0) cvm::fatal_error ("Error: providing an invalid number of points, "+ cvm::to_str (nx_i[i])+".\n"); nxc[i] = nt; nt *= nx[i]; } data.reserve (nt); data.assign (nt, t); } /// \brief Allocate data (allow initialization also after construction) void create() { create (this->nx, T(), this->mult); } /// \brief Reset data (in case the grid is being reused) void reset (T const &t = T()) { data.assign (nt, t); } /// Default constructor colvar_grid() : has_data (false) { save_delimiters = false; nd = nt = 0; } /// Destructor virtual ~colvar_grid() {} /// \brief "Almost copy-constructor": only copies configuration /// parameters from another grid, but doesn't reallocate stuff; /// create() must be called after that; colvar_grid (colvar_grid const &g) : has_data (false), nd (g.nd), + nx (g.nx), mult (g.mult), cv (g.cv), lower_boundaries (g.lower_boundaries), upper_boundaries (g.upper_boundaries), hard_lower_boundaries (g.hard_lower_boundaries), hard_upper_boundaries (g.hard_upper_boundaries), periodic (g.periodic), widths (g.widths), actual_value (g.actual_value), data() { save_delimiters = false; } /// \brief Constructor from explicit grid sizes \param nx_i Number /// of grid points along each dimension \param t Initial value for /// the function at each point (optional) \param mult_i Multiplicity /// of each value colvar_grid (std::vector const &nx_i, T const &t = T(), size_t const &mult_i = 1) : has_data (false) { save_delimiters = false; this->create (nx_i, t, mult_i); } /// \brief Constructor from a vector of colvars colvar_grid (std::vector const &colvars, T const &t = T(), size_t const &mult_i = 1, bool margin = false) : cv (colvars), has_data (false) { save_delimiters = false; std::vector nx_i; if (cvm::debug()) cvm::log ("Allocating a grid for "+cvm::to_str (colvars.size())+ " collective variables.\n"); for (size_t i = 0; i < cv.size(); i++) { if (cv[i]->type() != colvarvalue::type_scalar) { cvm::fatal_error ("Colvar grids can only be automatically " "constructed for scalar variables. " "ABF and histogram can not be used; metadynamics " "can be used with useGrids disabled.\n"); } if (cv[i]->width <= 0.0) { cvm::fatal_error ("Tried to initialize a grid on a " "variable with negative or zero width.\n"); } if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) { cvm::fatal_error ("Tried to initialize a grid on a " "variable with undefined boundaries.\n"); } widths.push_back (cv[i]->width); hard_lower_boundaries.push_back (cv[i]->hard_lower_boundary); hard_upper_boundaries.push_back (cv[i]->hard_upper_boundary); periodic.push_back (cv[i]->periodic_boundaries()); // By default, get reported colvar value (for extended Lagrangian colvars) actual_value.push_back (false); // except if a colvar is specified twice in a row // then the first instance is the actual value if (i > 0 && cv[i-1] == cv[i]) { actual_value[i-1] = true; } if (margin) { if (periodic[i]) { // Shift the grid by half the bin width (values at edges instead of center of bins) lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]); upper_boundaries.push_back (cv[i]->upper_boundary.real_value - 0.5 * widths[i]); } else { // Make this grid larger by one bin width lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]); upper_boundaries.push_back (cv[i]->upper_boundary.real_value + 0.5 * widths[i]); } } else { lower_boundaries.push_back (cv[i]->lower_boundary); upper_boundaries.push_back (cv[i]->upper_boundary); } { cvm::real nbins = ( upper_boundaries[i].real_value - lower_boundaries[i].real_value ) / widths[i]; int nbins_round = (int)(nbins+0.5); if (std::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) { cvm::log ("Warning: grid interval ("+ cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+ cvm::to_str (upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+ ") is not commensurate to its bin width ("+ cvm::to_str (widths[i], cvm::cv_width, cvm::cv_prec)+").\n"); upper_boundaries[i].real_value = lower_boundaries[i].real_value + (nbins_round * widths[i]); } if (cvm::debug()) cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+ " for the colvar no. "+cvm::to_str (i+1)+".\n"); nx_i.push_back (nbins_round); } } create (nx_i, t, mult_i); } /// Wrap an index vector around periodic boundary conditions /// also checks validity of non-periodic indices inline void wrap (std::vector & ix) const { for (size_t i = 0; i < nd; i++) { if (periodic[i]) { ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result } else { if (ix[i] < 0 || ix[i] >= nx[i]) cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): " + cvm::to_str (ix)); } } } /// \brief Report the bin corresponding to the current value of variable i inline int current_bin_scalar(int const i) const { return value_to_bin_scalar (actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); } /// \brief Use the lower boundary and the width to report which bin /// the provided value is in inline int value_to_bin_scalar (colvarvalue const &value, const int i) const { return (int) std::floor ( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); } /// \brief Same as the standard version, but uses another grid definition inline int value_to_bin_scalar (colvarvalue const &value, colvarvalue const &new_offset, cvm::real const &new_width) const { return (int) std::floor ( (value.real_value - new_offset.real_value) / new_width ); } /// \brief Use the two boundaries and the width to report the /// central value corresponding to a bin index inline colvarvalue bin_to_value_scalar (int const &i_bin, int const i) const { return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin); } /// \brief Same as the standard version, but uses different parameters inline colvarvalue bin_to_value_scalar (int const &i_bin, colvarvalue const &new_offset, cvm::real const &new_width) const { return new_offset.real_value + new_width * (0.5 + i_bin); } /// Set the value at the point with index ix inline void set_value (std::vector const &ix, T const &t, size_t const &imult = 0) { data[this->address (ix)+imult] = t; has_data = true; } /// \brief Get the binned value indexed by ix, or the first of them /// if the multiplicity is larger than 1 inline T const & value (std::vector const &ix, size_t const &imult = 0) const { return data[this->address (ix) + imult]; } /// \brief Add a constant to all elements (fast loop) inline void add_constant (T const &t) { for (size_t i = 0; i < nt; i++) data[i] += t; has_data = true; } /// \brief Multiply all elements by a scalar constant (fast loop) inline void multiply_constant (cvm::real const &a) { for (size_t i = 0; i < nt; i++) data[i] *= a; } /// \brief Get the bin indices corresponding to the provided values of /// the colvars inline std::vector const get_colvars_index (std::vector const &values) const { std::vector index = new_index(); for (size_t i = 0; i < nd; i++) { index[i] = value_to_bin_scalar (values[i], i); } return index; } /// \brief Get the bin indices corresponding to the current values /// of the colvars inline std::vector const get_colvars_index() const { std::vector index = new_index(); for (size_t i = 0; i < nd; i++) { index[i] = current_bin_scalar (i); } return index; } /// \brief Get the minimal distance (in number of bins) from the /// boundaries; a negative number is returned if the given point is /// off-grid inline cvm::real bin_distance_from_boundaries (std::vector const &values, bool skip_hard_boundaries = false) { cvm::real minimum = 1.0E+16; for (size_t i = 0; i < nd; i++) { if (periodic[i]) continue; cvm::real dl = std::sqrt (cv[i]->dist2 (values[i], lower_boundaries[i])) / widths[i]; cvm::real du = std::sqrt (cv[i]->dist2 (values[i], upper_boundaries[i])) / widths[i]; if (values[i].real_value < lower_boundaries[i]) dl *= -1.0; if (values[i].real_value > upper_boundaries[i]) du *= -1.0; if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum)) minimum = dl; if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum)) minimum = du; } return minimum; } /// \brief Add data from another grid of the same type /// /// Note: this function maps other_grid inside this one regardless /// of whether it fits or not. void map_grid (colvar_grid const &other_grid) { if (other_grid.multiplicity() != this->multiplicity()) cvm::fatal_error ("Error: trying to merge two grids with values of " "different multiplicity.\n"); std::vector const &gb = this->lower_boundaries; std::vector const &gw = this->widths; std::vector const &ogb = other_grid.lower_boundaries; std::vector const &ogw = other_grid.widths; std::vector ix = this->new_index(); std::vector oix = other_grid.new_index(); if (cvm::debug()) cvm::log ("Remapping grid...\n"); for ( ; this->index_ok (ix); this->incr (ix)) { for (size_t i = 0; i < nd; i++) { oix[i] = value_to_bin_scalar (bin_to_value_scalar (ix[i], gb[i], gw[i]), ogb[i], ogw[i]); } if (! other_grid.index_ok (oix)) { continue; } for (size_t im = 0; im < mult; im++) { this->set_value (ix, other_grid.value (oix, im), im); } } has_data = true; if (cvm::debug()) cvm::log ("Remapping done.\n"); } /// \brief Add data from another grid of the same type, AND /// identical definition (boundaries, widths) void add_grid (colvar_grid const &other_grid, cvm::real scale_factor = 1.0) { if (other_grid.multiplicity() != this->multiplicity()) cvm::fatal_error ("Error: trying to sum togetehr two grids with values of " "different multiplicity.\n"); if (scale_factor != 1.0) for (size_t i = 0; i < data.size(); i++) { data[i] += scale_factor * other_grid.data[i]; } else // skip multiplication if possible for (size_t i = 0; i < data.size(); i++) { data[i] += other_grid.data[i]; } has_data = true; } /// \brief Return the value suitable for output purposes (so that it /// may be rescaled or manipulated without changing it permanently) virtual inline T value_output (std::vector const &ix, size_t const &imult = 0) { return value (ix, imult); } /// \brief Get the value from a formatted output and transform it /// into the internal representation (the two may be different, /// e.g. when using colvar_grid_count) virtual inline void value_input (std::vector const &ix, T const &t, size_t const &imult = 0, bool add = false) { if ( add ) data[address (ix) + imult] += t; else data[address (ix) + imult] = t; has_data = true; } // /// Get the pointer to the binned value indexed by ix // inline T const *value_p (std::vector const &ix) // { // return &(data[address (ix)]); // } /// \brief Get the index corresponding to the "first" bin, to be /// used as the initial value for an index in looping inline std::vector const new_index() const { return std::vector (nd, 0); } /// \brief Check that the index is within range in each of the /// dimensions inline bool index_ok (std::vector const &ix) const { for (size_t i = 0; i < nd; i++) { if ( (ix[i] < 0) || (ix[i] >= int (nx[i])) ) return false; } return true; } /// \brief Increment the index, in a way that will make it loop over /// the whole nd-dimensional array inline void incr (std::vector &ix) const { for (int i = ix.size()-1; i >= 0; i--) { ix[i]++; if (ix[i] >= nx[i]) { if (i > 0) { ix[i] = 0; continue; } else { // this is the last iteration, a non-valid index is being // set for the outer index, which will be caught by // index_ok() ix[0] = nx[0]; return; } } else { return; } } } /// \brief Write the grid parameters (number of colvars, boundaries, width and number of points) std::ostream & write_params (std::ostream &os) { os << "grid_parameters {\n n_colvars " << nd << "\n"; os << " lower_boundaries "; for (size_t i = 0; i < nd; i++) os << " " << lower_boundaries[i]; os << "\n"; os << " upper_boundaries "; for (size_t i = 0; i < nd; i++) os << " " << upper_boundaries[i]; os << "\n"; os << " widths "; for (size_t i = 0; i < nd; i++) os << " " << widths[i]; os << "\n"; os << " sizes "; for (size_t i = 0; i < nd; i++) os << " " << nx[i]; os << "\n"; os << "}\n"; return os; } bool parse_params (std::string const &conf) { std::vector old_nx = nx; std::vector old_lb = lower_boundaries; { size_t nd_in = 0; colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent); if (nd_in != nd) cvm::fatal_error ("Error: trying to read data for a grid " "that contains a different number of colvars ("+ cvm::to_str (nd_in)+") than the grid defined " "in the configuration file ("+cvm::to_str (nd)+ ").\n"); } colvarparse::get_keyval (conf, "lower_boundaries", lower_boundaries, lower_boundaries, colvarparse::parse_silent); colvarparse::get_keyval (conf, "upper_boundaries", upper_boundaries, upper_boundaries, colvarparse::parse_silent); colvarparse::get_keyval (conf, "widths", widths, widths, colvarparse::parse_silent); colvarparse::get_keyval (conf, "sizes", nx, nx, colvarparse::parse_silent); bool new_params = false; for (size_t i = 0; i < nd; i++) { if ( (old_nx[i] != nx[i]) || (std::sqrt (cv[i]->dist2 (old_lb[i], lower_boundaries[i])) > 1.0E-10) ) { new_params = true; } } // reallocate the array in case the grid params have just changed if (new_params) { data.resize (0); this->create (nx, T(), mult); } return true; } /// \brief Check that the grid information inside (boundaries, /// widths, ...) is consistent with the current setting of the /// colvars void check_consistency() { for (size_t i = 0; i < nd; i++) { if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary, lower_boundaries[i])) > 1.0E-10) || (std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary, upper_boundaries[i])) > 1.0E-10) || (std::sqrt (cv[i]->dist2 (cv[i]->width, widths[i])) > 1.0E-10) ) { cvm::fatal_error ("Error: restart information for a grid is " "inconsistent with that of its colvars.\n"); } } } /// \brief Check that the grid information inside (boundaries, /// widths, ...) is consistent with the one of another grid void check_consistency (colvar_grid const &other_grid) { for (size_t i = 0; i < nd; i++) { // we skip dist2(), because periodicities and the like should // matter: boundaries should be EXACTLY the same (otherwise, // map_grid() should be used) if ( (std::fabs (other_grid.lower_boundaries[i] - lower_boundaries[i]) > 1.0E-10) || (std::fabs (other_grid.upper_boundaries[i] - upper_boundaries[i]) > 1.0E-10) || (std::fabs (other_grid.widths[i] - widths[i]) > 1.0E-10) || (data.size() != other_grid.data.size()) ) { cvm::fatal_error ("Error: inconsistency between " "two grids that are supposed to be equal, " "aside from the data stored.\n"); } } } /// \brief Write the grid data without labels, as they are /// represented in memory /// \param buf_size Number of values per line std::ostream & write_raw (std::ostream &os, size_t const buf_size = 3) { std::streamsize const w = os.width(); std::streamsize const p = os.precision(); std::vector ix = new_index(); size_t count = 0; for ( ; index_ok (ix); incr (ix)) { for (size_t imult = 0; imult < mult; imult++) { os << " " << std::setw (w) << std::setprecision (p) << value_output (ix, imult); if (((++count) % buf_size) == 0) os << "\n"; } } // write a final newline only if buffer is not empty if ((count % buf_size) != 0) os << "\n"; return os; } /// \brief Read data written by colvar_grid::write_raw() std::istream & read_raw (std::istream &is) { size_t const start_pos = is.tellg(); for (std::vector ix = new_index(); index_ok (ix); incr (ix)) { for (size_t imult = 0; imult < mult; imult++) { T new_value; if (is >> new_value) { value_input (ix, new_value, imult); } else { is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } } } has_data = true; return is; } /// \brief To be called after colvar_grid::read_raw() returns an error void read_raw_error() { cvm::fatal_error ("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n"); } /// \brief Write the grid in a format which is both human readable /// and suitable for visualization e.g. with gnuplot void write_multicol (std::ostream &os) { std::streamsize const w = os.width(); std::streamsize const p = os.precision(); // Data in the header: nColvars, then for each // xiMin, dXi, nPoints, periodic os << std::setw (2) << "# " << nd << "\n"; for (size_t i = 0; i < nd; i++) { os << "# " << std::setw (10) << lower_boundaries[i] << std::setw (10) << widths[i] << std::setw (10) << nx[i] << " " << periodic[i] << "\n"; } for (std::vector ix = new_index(); index_ok (ix); incr (ix) ) { if (ix.back() == 0) { // if the last index is 0, add a new line to mark the new record os << "\n"; } for (size_t i = 0; i < nd; i++) { os << " " << std::setw (w) << std::setprecision (p) << bin_to_value_scalar (ix[i], i); } os << " "; for (size_t imult = 0; imult < mult; imult++) { os << " " << std::setw (w) << std::setprecision (p) << value_output (ix, imult); } os << "\n"; } } /// \brief Read a grid written by colvar_grid::write_multicol() /// Adding data if add is true, replacing if false std::istream & read_multicol (std::istream &is, bool add = false) { // Data in the header: nColvars, then for each // xiMin, dXi, nPoints, periodic std::string hash; cvm::real lower, width, x; size_t n, periodic; bool remap; std::vector new_value; std::vector nx_read; std::vector bin; if ( cv.size() != nd ) { cvm::fatal_error ("Cannot read grid file: missing reference to colvars."); } if ( !(is >> hash) || (hash != "#") ) { cvm::fatal_error ("Error reading grid at position "+ cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n"); } is >> n; if ( n != nd ) { cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n"); } nx_read.resize (n); bin.resize (n); new_value.resize (mult); if (this->has_parent_data && add) { new_data.resize (data.size()); } remap = false; for (size_t i = 0; i < nd; i++ ) { if ( !(is >> hash) || (hash != "#") ) { cvm::fatal_error ("Error reading grid at position "+ cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n"); } is >> lower >> width >> nx_read[i] >> periodic; if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) || (std::fabs (width - widths[i] ) > 1.0e-10) || (nx_read[i] != nx[i]) ) { cvm::log ("Warning: reading from different grid definition (colvar " + cvm::to_str (i+1) + "); remapping data on new grid.\n"); remap = true; } } if ( remap ) { // re-grid data while (is.good()) { bool end_of_file = false; for (size_t i = 0; i < nd; i++ ) { if ( !(is >> x) ) end_of_file = true; bin[i] = value_to_bin_scalar (x, i); } if (end_of_file) break; for (size_t imult = 0; imult < mult; imult++) { is >> new_value[imult]; } if ( index_ok(bin) ) { for (size_t imult = 0; imult < mult; imult++) { value_input (bin, new_value[imult], imult, add); } } } } else { // do not re-grid the data but assume the same grid is used for (std::vector ix = new_index(); index_ok (ix); incr (ix) ) { for (size_t i = 0; i < nd; i++ ) { is >> x; } for (size_t imult = 0; imult < mult; imult++) { is >> new_value[imult]; value_input (ix, new_value[imult], imult, add); } } } has_data = true; return is; } }; /// \brief Colvar_grid derived class to hold counters in discrete /// n-dim colvar space class colvar_grid_count : public colvar_grid { public: /// Default constructor colvar_grid_count(); /// Destructor virtual inline ~colvar_grid_count() {} /// Constructor colvar_grid_count (std::vector const &nx_i, size_t const &def_count = 0); /// Constructor from a vector of colvars colvar_grid_count (std::vector &colvars, size_t const &def_count = 0); /// Increment the counter at given position inline void incr_count (std::vector const &ix) { ++(data[this->address (ix)]); } /// \brief Get the binned count indexed by ix from the newly read data inline size_t const & new_count (std::vector const &ix, size_t const &imult = 0) { return new_data[address (ix) + imult]; } /// \brief Read the grid from a restart std::istream & read_restart (std::istream &is); /// \brief Write the grid to a restart std::ostream & write_restart (std::ostream &os); /// \brief Get the value from a formatted output and transform it /// into the internal representation (it may have been rescaled or /// manipulated) virtual inline void value_input (std::vector const &ix, size_t const &t, size_t const &imult = 0, bool add = false) { if (add) { data[address (ix)] += t; if (this->has_parent_data) { // save newly read data for inputting parent grid new_data[address (ix)] = t; } } else { data[address (ix)] = t; } has_data = true; } }; /// Class for accumulating a scalar function on a grid class colvar_grid_scalar : public colvar_grid { public: /// \brief Provide the associated sample count by which each binned value /// should be divided colvar_grid_count *samples; /// Default constructor colvar_grid_scalar(); /// Copy constructor (needed because of the grad pointer) colvar_grid_scalar (colvar_grid_scalar const &g); /// Destructor ~colvar_grid_scalar(); /// Constructor from specific sizes arrays colvar_grid_scalar (std::vector const &nx_i); /// Constructor from a vector of colvars colvar_grid_scalar (std::vector &colvars, bool margin = 0); /// Accumulate the value inline void acc_value (std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0) { // only legal value of imult here is 0 data[address (ix)] += new_value; if (samples) samples->incr_count (ix); has_data = true; } /// Return the gradient of the scalar field from finite differences inline const cvm::real * gradient_finite_diff ( const std::vector &ix0 ) { cvm::real A0, A1; std::vector ix; if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only."); for (int n = 0; n < nd; n++) { ix = ix0; A0 = data[address (ix)]; ix[n]++; wrap (ix); A1 = data[address (ix)]; ix[1-n]++; wrap (ix); A1 += data[address (ix)]; ix[n]--; wrap (ix); A0 += data[address (ix)]; grad[n] = 0.5 * (A1 - A0) / widths[n]; } return grad; } /// \brief Return the value of the function at ix divided by its /// number of samples (if the count grid is defined) virtual cvm::real value_output (std::vector const &ix, size_t const &imult = 0) { if (imult > 0) cvm::fatal_error ("Error: trying to access a component " "larger than 1 in a scalar data grid.\n"); if (samples) return (samples->value (ix) > 0) ? (data[address (ix)] / cvm::real (samples->value (ix))) : 0.0; else return data[address (ix)]; } /// \brief Get the value from a formatted output and transform it /// into the internal representation (it may have been rescaled or /// manipulated) virtual void value_input (std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0, bool add = false) { if (imult > 0) cvm::fatal_error ("Error: trying to access a component " "larger than 1 in a scalar data grid.\n"); if (add) { if (samples) data[address (ix)] += new_value * samples->new_count (ix); else data[address (ix)] += new_value; } else { if (samples) data[address (ix)] = new_value * samples->value (ix); else data[address (ix)] = new_value; } has_data = true; } /// \brief Read the grid from a restart std::istream & read_restart (std::istream &is); /// \brief Write the grid to a restart std::ostream & write_restart (std::ostream &os); /// \brief Return the highest value inline cvm::real maximum_value() { cvm::real max = data[0]; for (size_t i = 0; i < nt; i++) { if (data[i] > max) max = data[i]; } return max; } /// \brief Return the lowest value inline cvm::real minimum_value() { cvm::real min = data[0]; for (size_t i = 0; i < nt; i++) { if (data[i] < min) min = data[i]; } return min; } private: // gradient cvm::real * grad; }; /// Class for accumulating the gradient of a scalar function on a grid class colvar_grid_gradient : public colvar_grid { public: /// \brief Provide the sample count by which each binned value /// should be divided colvar_grid_count *samples; /// Default constructor colvar_grid_gradient(); /// Destructor virtual inline ~colvar_grid_gradient() {} /// Constructor from specific sizes arrays colvar_grid_gradient (std::vector const &nx_i); /// Constructor from a vector of colvars colvar_grid_gradient (std::vector &colvars); /// \brief Accumulate the gradient inline void acc_grad (std::vector const &ix, cvm::real const *grads) { for (size_t imult = 0; imult < mult; imult++) { data[address (ix) + imult] += grads[imult]; } if (samples) samples->incr_count (ix); } /// \brief Accumulate the gradient based on the force (i.e. sums the /// opposite of the force) inline void acc_force (std::vector const &ix, cvm::real const *forces) { for (size_t imult = 0; imult < mult; imult++) { data[address (ix) + imult] -= forces[imult]; } if (samples) samples->incr_count (ix); } /// \brief Return the value of the function at ix divided by its /// number of samples (if the count grid is defined) virtual inline cvm::real value_output (std::vector const &ix, size_t const &imult = 0) { if (samples) return (samples->value (ix) > 0) ? (data[address (ix) + imult] / cvm::real (samples->value (ix))) : 0.0; else return data[address (ix) + imult]; } /// \brief Get the value from a formatted output and transform it /// into the internal representation (it may have been rescaled or /// manipulated) virtual inline void value_input (std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0, bool add = false) { if (add) { if (samples) data[address (ix) + imult] += new_value * samples->new_count (ix); else data[address (ix) + imult] += new_value; } else { if (samples) data[address (ix) + imult] = new_value * samples->value (ix); else data[address (ix) + imult] = new_value; } has_data = true; } /// \brief Read the grid from a restart std::istream & read_restart (std::istream &is); /// \brief Write the grid to a restart std::ostream & write_restart (std::ostream &os); /// Compute and return average value for a 1D gradient grid inline cvm::real average() { size_t n = 0; if (nd != 1 || nx[0] == 0) { return 0.0; } cvm::real sum = 0.0; std::vector ix = new_index(); if (samples) { for ( ; index_ok (ix); incr (ix)) { if ( (n = samples->value (ix)) ) sum += value (ix) / n; } } else { for ( ; index_ok (ix); incr (ix)) { sum += value (ix); } } return (sum / cvm::real (nx[0])); } /// \brief If the grid is 1-dimensional, integrate it and write the /// integral to a file void write_1D_integral (std::ostream &os); }; #endif diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 00744740c..67ef6cadd 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1,1369 +1,1369 @@ #include "colvarmodule.h" #include "colvarparse.h" #include "colvarproxy.h" #include "colvar.h" #include "colvarbias.h" #include "colvarbias_meta.h" #include "colvarbias_abf.h" colvarmodule::colvarmodule (char const *config_filename, colvarproxy *proxy_in) { // pointer to the proxy object if (proxy == NULL) { proxy = proxy_in; parse = new colvarparse(); } else { cvm::fatal_error ("Error: trying to allocate twice the collective " "variable module.\n"); } cvm::log (cvm::line_marker); cvm::log ("Initializing the collective variables module, version "+ cvm::to_str(COLVARS_VERSION)+".\n"); // "it_restart" will be set by the input state file, if any; // "it" should be updated by the proxy it = it_restart = 0; it_restart_from_state_file = true; // open the configfile config_s.open (config_filename); if (!config_s) cvm::fatal_error ("Error: in opening configuration file \""+ std::string (config_filename)+"\".\n"); // read the config file into a string std::string conf = ""; { std::string line; while (colvarparse::getline_nocomments (config_s, line)) conf.append (line+"\n"); // don't need the stream any more config_s.close(); } parse->get_keyval (conf, "analysis", b_analysis, false); - parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03, + parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07, colvarparse::parse_silent); parse->get_keyval (conf, "eigenvalueCrossingThreshold", - colvarmodule::rotation::crossing_threshold, 1.0e-04, + colvarmodule::rotation::crossing_threshold, 1.0e-02, colvarparse::parse_silent); parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100); parse->get_keyval (conf, "colvarsRestartFrequency", restart_out_freq, proxy->restart_frequency()); // by default overwrite the existing trajectory file parse->get_keyval (conf, "colvarsTrajAppend", cv_traj_append, false); // input restart file restart_in_name = proxy->input_prefix().size() ? std::string (proxy->input_prefix()+".colvars.state") : std::string ("") ; // output restart file restart_out_name = proxy->restart_output_prefix().size() ? std::string (proxy->restart_output_prefix()+".colvars.state") : std::string (""); if (restart_out_name.size()) cvm::log ("The restart output state file will be \""+restart_out_name+"\".\n"); output_prefix = proxy->output_prefix(); cvm::log ("The final output state file will be \""+ (output_prefix.size() ? std::string (output_prefix+".colvars.state") : std::string ("colvars.state"))+"\".\n"); cv_traj_name = (output_prefix.size() ? std::string (output_prefix+".colvars.traj") : std::string ("colvars.traj")); cvm::log ("The trajectory file will be \""+ cv_traj_name+"\".\n"); if (cv_traj_freq) { // open trajectory file if (cv_traj_append) { cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+ "\".\n"); cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); } else { proxy->backup_file (cv_traj_name.c_str()); cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); } cv_traj_os.setf (std::ios::scientific, std::ios::floatfield); } // parse the options for collective variables init_colvars (conf); // parse the options for biases init_biases (conf); // done with the parsing, check that all keywords are valid parse->check_keywords (conf, "colvarmodule"); cvm::log (cvm::line_marker); // read the restart configuration, if available if (restart_in_name.size()) { // read the restart file std::ifstream input_is (restart_in_name.c_str()); if (!input_is.good()) fatal_error ("Error: in opening restart file \""+ std::string (restart_in_name)+"\".\n"); else { cvm::log ("Restarting from file \""+restart_in_name+"\".\n"); read_restart (input_is); cvm::log (cvm::line_marker); } } // check if it is possible to save output configuration if ((!output_prefix.size()) && (!restart_out_name.size())) { cvm::fatal_error ("Error: neither the final output state file or " "the output restart file could be defined, exiting.\n"); } cvm::log ("Collective variables module initialized.\n"); cvm::log (cvm::line_marker); } std::istream & colvarmodule::read_restart (std::istream &is) { { // read global restart information std::string restart_conf; if (is >> colvarparse::read_block ("configuration", restart_conf)) { if (it_restart_from_state_file) { parse->get_keyval (restart_conf, "step", it_restart, (size_t) 0, colvarparse::parse_silent); it = it_restart; } } is.clear(); } // colvars restart cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ( !((*cvi)->read_restart (is)) ) cvm::fatal_error ("Error: in reading restart configuration for collective variable \""+ (*cvi)->name+"\".\n"); } // biases restart for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if (!((*bi)->read_restart (is))) fatal_error ("Error: in reading restart configuration for bias \""+ (*bi)->name+"\".\n"); } cvm::decrease_depth(); return is; } void colvarmodule::init_colvars (std::string const &conf) { if (cvm::debug()) cvm::log ("Initializing the collective variables.\n"); std::string colvar_conf = ""; size_t pos = 0; while (parse->key_lookup (conf, "colvar", colvar_conf, pos)) { if (colvar_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); colvars.push_back (new colvar (colvar_conf)); (colvars.back())->check_keywords (colvar_conf, "colvar"); cvm::decrease_depth(); } else { cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n"); } colvar_conf = ""; } if (!colvars.size()) cvm::fatal_error ("Error: no collective variables defined.\n"); if (colvars.size()) cvm::log (cvm::line_marker); cvm::log ("Collective variables initialized, "+ cvm::to_str (colvars.size())+ " in total.\n"); } void colvarmodule::init_biases (std::string const &conf) { if (cvm::debug()) cvm::log ("Initializing the collective variables biases.\n"); { /// initialize ABF instances std::string abf_conf = ""; size_t abf_pos = 0; while (parse->key_lookup (conf, "abf", abf_conf, abf_pos)) { if (abf_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_abf (abf_conf, "abf")); (biases.back())->check_keywords (abf_conf, "abf"); cvm::decrease_depth(); n_abf_biases++; } else { cvm::log ("Warning: \"abf\" keyword found without configuration.\n"); } abf_conf = ""; } } { /// initialize harmonic restraints std::string harm_conf = ""; size_t harm_pos = 0; while (parse->key_lookup (conf, "harmonic", harm_conf, harm_pos)) { if (harm_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_harmonic (harm_conf, "harmonic")); (biases.back())->check_keywords (harm_conf, "harmonic"); cvm::decrease_depth(); n_harm_biases++; } else { cvm::log ("Warning: \"harmonic\" keyword found without configuration.\n"); } harm_conf = ""; } } { /// initialize histograms std::string histo_conf = ""; size_t histo_pos = 0; while (parse->key_lookup (conf, "histogram", histo_conf, histo_pos)) { if (histo_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_histogram (histo_conf, "histogram")); (biases.back())->check_keywords (histo_conf, "histogram"); cvm::decrease_depth(); n_histo_biases++; } else { cvm::log ("Warning: \"histogram\" keyword found without configuration.\n"); } histo_conf = ""; } } { /// initialize metadynamics instances std::string meta_conf = ""; size_t meta_pos = 0; while (parse->key_lookup (conf, "metadynamics", meta_conf, meta_pos)) { if (meta_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_meta (meta_conf, "metadynamics")); (biases.back())->check_keywords (meta_conf, "metadynamics"); cvm::decrease_depth(); n_meta_biases++; } else { cvm::log ("Warning: \"metadynamics\" keyword found without configuration.\n"); } meta_conf = ""; } } if (biases.size()) cvm::log (cvm::line_marker); cvm::log ("Collective variables biases initialized, "+ cvm::to_str (biases.size())+" in total.\n"); } void colvarmodule::change_configuration( std::string const &name, std::string const &conf) { cvm::increase_depth(); int found = 0; for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if ( (*bi)->name == name ) { ++found; (*bi)->change_configuration(conf); } } if ( found < 1 ) cvm::fatal_error ("Error: bias not found"); if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name"); cvm::decrease_depth(); } cvm::real colvarmodule::energy_difference( std::string const &name, std::string const &conf) { cvm::increase_depth(); cvm::real energy_diff = 0.; int found = 0; for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if ( (*bi)->name == name ) { ++found; energy_diff = (*bi)->energy_difference(conf); } } if ( found < 1 ) cvm::fatal_error ("Error: bias not found"); if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name"); cvm::decrease_depth(); return energy_diff; } void colvarmodule::calc() { cvm::real total_bias_energy = 0.0; cvm::real total_colvar_energy = 0.0; if (cvm::debug()) { cvm::log (cvm::line_marker); cvm::log ("Collective variables module, step no. "+ cvm::to_str (cvm::step_absolute())+"\n"); } // calculate collective variables and their gradients if (cvm::debug()) cvm::log ("Calculating collective variables.\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->calc(); } cvm::decrease_depth(); // update the biases and communicate their forces to the collective // variables if (cvm::debug() && biases.size()) cvm::log ("Updating collective variable biases.\n"); cvm::increase_depth(); for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { total_bias_energy += (*bi)->update(); } cvm::decrease_depth(); // sum the forces from all biases for each collective variable if (cvm::debug() && biases.size()) cvm::log ("Collecting forces from all biases.\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->reset_bias_force(); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->communicate_forces(); } cvm::decrease_depth(); if (cvm::b_analysis) { // perform runtime analysis of colvars and biases if (cvm::debug() && biases.size()) cvm::log ("Perform runtime analyses.\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->analyse(); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->analyse(); } cvm::decrease_depth(); } // sum up the forces for each colvar and integrate any internal // equation of motion if (cvm::debug()) cvm::log ("Updating the internal degrees of freedom " "of colvars (if they have any).\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { total_colvar_energy += (*cvi)->update(); } cvm::decrease_depth(); proxy->add_energy (total_bias_energy + total_colvar_energy); // make collective variables communicate their forces to their // coupled degrees of freedom (i.e. atoms) if (cvm::debug()) cvm::log ("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ((*cvi)->tasks[colvar::task_gradients]) (*cvi)->communicate_forces(); } cvm::decrease_depth(); // write restart file, if needed if (restart_out_freq && !cvm::b_analysis) { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log ("Writing the state file \""+ restart_out_name+"\".\n"); proxy->backup_file (restart_out_name.c_str()); restart_out_os.open (restart_out_name.c_str()); restart_out_os.setf (std::ios::scientific, std::ios::floatfield); if (!write_restart (restart_out_os)) cvm::fatal_error ("Error: in writing restart file.\n"); restart_out_os.close(); } } // write trajectory file, if needed if (cv_traj_freq) { if (cvm::debug()) cvm::log ("Writing trajectory file.\n"); // (re)open trajectory file if (!cv_traj_os.good()) { if (cv_traj_append) { cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+ "\".\n"); cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); } else { cvm::log ("Overwriting colvar trajectory file \""+cv_traj_name+ "\".\n"); proxy->backup_file (cv_traj_name.c_str()); cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); } cv_traj_os.setf (std::ios::scientific, std::ios::floatfield); } // write labels in the traj file every 1000 lines and at first ts cvm::increase_depth(); if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { cv_traj_os << "# " << cvm::wrap_string ("step", cvm::it_width-2) << " "; if (cvm::debug()) cv_traj_os.flush(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj_label (cv_traj_os); } cv_traj_os << "\n"; if (cvm::debug()) cv_traj_os.flush(); } cvm::decrease_depth(); // write collective variable values to the traj file cvm::increase_depth(); if ((cvm::step_absolute() % cv_traj_freq) == 0) { cv_traj_os << std::setw (cvm::it_width) << it << " "; for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj (cv_traj_os); } cv_traj_os << "\n"; if (cvm::debug()) cv_traj_os.flush(); } cvm::decrease_depth(); if (restart_out_freq) { // flush the trajectory file if we are at the restart frequency if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log ("Synchronizing (emptying the buffer of) trajectory file \""+ cv_traj_name+"\".\n"); cv_traj_os.flush(); } } } // end if (cv_traj_freq) } void colvarmodule::analyze() { if (cvm::debug()) { cvm::log ("colvarmodule::analyze(), step = "+cvm::to_str (it)+".\n"); } if (cvm::step_relative() == 0) cvm::log ("Performing analysis.\n"); // perform colvar-specific analysis for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { cvm::increase_depth(); (*cvi)->analyse(); cvm::decrease_depth(); } // perform bias-specific analysis for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { cvm::increase_depth(); (*bi)->analyse(); cvm::decrease_depth(); } } colvarmodule::~colvarmodule() { for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { delete *cvi; } colvars.clear(); for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { delete *bi; } biases.clear(); if (cv_traj_os.good()) { cv_traj_os.close(); } delete parse; proxy = NULL; } void colvarmodule::write_output_files() { // if this is a simulation run (i.e. not a postprocessing), output data // must be written to be able to restart the simulation std::string const out_name = (output_prefix.size() ? std::string (output_prefix+".colvars.state") : std::string ("colvars.state")); cvm::log ("Saving collective variables state to \""+out_name+"\".\n"); proxy->backup_file (out_name.c_str()); std::ofstream out (out_name.c_str()); out.setf (std::ios::scientific, std::ios::floatfield); this->write_restart (out); out.close(); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_output_files(); } cvm::decrease_depth(); // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); } bool colvarmodule::read_traj (char const *traj_filename, size_t traj_read_begin, size_t traj_read_end) { cvm::log ("Opening trajectory file \""+ std::string (traj_filename)+"\".\n"); std::ifstream traj_is (traj_filename); while (true) { while (true) { std::string line (""); do { if (!colvarparse::getline_nocomments (traj_is, line)) { cvm::log ("End of file \""+std::string (traj_filename)+ "\" reached, or corrupted file.\n"); traj_is.close(); return false; } } while (line.find_first_not_of (colvarparse::white_space) == std::string::npos); std::istringstream is (line); if (!(is >> it)) return false; if ( (it < traj_read_begin) ) { if ((it % 1000) == 0) std::cerr << "Skipping trajectory step " << it << " \r"; continue; } else { if ((it % 1000) == 0) std::cerr << "Reading from trajectory, step = " << it << " \r"; if ( (traj_read_end > traj_read_begin) && (it > traj_read_end) ) { std::cerr << "\n"; cvm::log ("Reached the end of the trajectory, " "read_end = "+cvm::to_str (traj_read_end)+"\n"); return false; } for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->read_traj (is)) { cvm::log ("Error: in reading colvar \""+(*cvi)->name+ "\" from trajectory file \""+ std::string (traj_filename)+"\".\n"); return false; } } break; } } } return true; } std::ostream & colvarmodule::write_restart (std::ostream &os) { os << "configuration {\n" << " step " << std::setw (it_width) << it << "\n" << " dt " << dt() << "\n" << "}\n\n"; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_restart (os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_restart (os); } cvm::decrease_depth(); return os; } void cvm::log (std::string const &message) { if (depth > 0) proxy->log ((std::string (2*depth, ' '))+message); else proxy->log (message); } void cvm::increase_depth() { depth++; } void cvm::decrease_depth() { if (depth) depth--; } void cvm::fatal_error (std::string const &message) { proxy->fatal_error (message); } void cvm::exit (std::string const &message) { proxy->exit (message); } // static pointers std::vector colvarmodule::colvars; std::vector colvarmodule::biases; size_t colvarmodule::n_abf_biases = 0; size_t colvarmodule::n_harm_biases = 0; size_t colvarmodule::n_histo_biases = 0; size_t colvarmodule::n_meta_biases = 0; colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03; size_t colvarmodule::it = 0; size_t colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::depth = 0; bool colvarmodule::b_analysis = false; cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04; // file name prefixes std::string colvarmodule::output_prefix = ""; std::string colvarmodule::input_prefix = ""; std::string colvarmodule::restart_in_name = ""; // i/o constants size_t const colvarmodule::it_width = 12; size_t const colvarmodule::cv_prec = 14; size_t const colvarmodule::cv_width = 21; size_t const colvarmodule::en_prec = 14; size_t const colvarmodule::en_width = 21; std::string const colvarmodule::line_marker = "----------------------------------------------------------------------\n"; std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v) { std::streamsize const w = os.width(); std::streamsize const p = os.precision(); os.width (2); os << "( "; os.width (w); os.precision (p); os << v.x << " , "; os.width (w); os.precision (p); os << v.y << " , "; os.width (w); os.precision (p); os << v.z << " )"; return os; } std::istream & operator >> (std::istream &is, colvarmodule::rvector &v) { size_t const start_pos = is.tellg(); char sep; if ( !(is >> sep) || !(sep == '(') || !(is >> v.x) || !(is >> sep) || !(sep == ',') || !(is >> v.y) || !(is >> sep) || !(sep == ',') || !(is >> v.z) || !(is >> sep) || !(sep == ')') ) { is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } return is; } std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q) { std::streamsize const w = os.width(); std::streamsize const p = os.precision(); os.width (2); os << "( "; os.width (w); os.precision (p); os << q.q0 << " , "; os.width (w); os.precision (p); os << q.q1 << " , "; os.width (w); os.precision (p); os << q.q2 << " , "; os.width (w); os.precision (p); os << q.q3 << " )"; return os; } std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q) { size_t const start_pos = is.tellg(); std::string euler (""); if ( (is >> euler) && (colvarparse::to_lower_cppstr (euler) == std::string ("euler")) ) { // parse the Euler angles char sep; cvm::real phi, theta, psi; if ( !(is >> sep) || !(sep == '(') || !(is >> phi) || !(is >> sep) || !(sep == ',') || !(is >> theta) || !(is >> sep) || !(sep == ',') || !(is >> psi) || !(is >> sep) || !(sep == ')') ) { is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } q = colvarmodule::quaternion (phi, theta, psi); } else { // parse the quaternion components is.seekg (start_pos, std::ios::beg); char sep; if ( !(is >> sep) || !(sep == '(') || !(is >> q.q0) || !(is >> sep) || !(sep == ',') || !(is >> q.q1) || !(is >> sep) || !(sep == ',') || !(is >> q.q2) || !(is >> sep) || !(sep == ',') || !(is >> q.q3) || !(is >> sep) || !(sep == ')') ) { is.clear(); is.seekg (start_pos, std::ios::beg); is.setstate (std::ios::failbit); return is; } } return is; } cvm::quaternion cvm::quaternion::position_derivative_inner (cvm::rvector const &pos, cvm::rvector const &vec) const { cvm::quaternion result (0.0, 0.0, 0.0, 0.0); result.q0 = 2.0 * pos.x * q0 * vec.x +2.0 * pos.y * q0 * vec.y +2.0 * pos.z * q0 * vec.z -2.0 * pos.y * q3 * vec.x +2.0 * pos.z * q2 * vec.x +2.0 * pos.x * q3 * vec.y -2.0 * pos.z * q1 * vec.y -2.0 * pos.x * q2 * vec.z +2.0 * pos.y * q1 * vec.z; result.q1 = +2.0 * pos.x * q1 * vec.x -2.0 * pos.y * q1 * vec.y -2.0 * pos.z * q1 * vec.z +2.0 * pos.y * q2 * vec.x +2.0 * pos.z * q3 * vec.x +2.0 * pos.x * q2 * vec.y -2.0 * pos.z * q0 * vec.y +2.0 * pos.x * q3 * vec.z +2.0 * pos.y * q0 * vec.z; result.q2 = -2.0 * pos.x * q2 * vec.x +2.0 * pos.y * q2 * vec.y -2.0 * pos.z * q2 * vec.z +2.0 * pos.y * q1 * vec.x +2.0 * pos.z * q0 * vec.x +2.0 * pos.x * q1 * vec.y +2.0 * pos.z * q3 * vec.y -2.0 * pos.x * q0 * vec.z +2.0 * pos.y * q3 * vec.z; result.q3 = -2.0 * pos.x * q3 * vec.x -2.0 * pos.y * q3 * vec.y +2.0 * pos.z * q3 * vec.z -2.0 * pos.y * q0 * vec.x +2.0 * pos.z * q1 * vec.x +2.0 * pos.x * q0 * vec.y +2.0 * pos.z * q2 * vec.y +2.0 * pos.x * q1 * vec.z +2.0 * pos.y * q2 * vec.z; return result; } // Calculate the optimal rotation between two groups, and implement it // as a quaternion. The method is the one documented in: Coutsias EA, // Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput // Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254 void colvarmodule::rotation::build_matrix (std::vector const &pos1, std::vector const &pos2, matrix2d &S) { cvm::rmatrix C; // build the correlation matrix C.reset(); for (size_t i = 0; i < pos1.size(); i++) { C.xx() += pos1[i].x * pos2[i].x; C.xy() += pos1[i].x * pos2[i].y; C.xz() += pos1[i].x * pos2[i].z; C.yx() += pos1[i].y * pos2[i].x; C.yy() += pos1[i].y * pos2[i].y; C.yz() += pos1[i].y * pos2[i].z; C.zx() += pos1[i].z * pos2[i].x; C.zy() += pos1[i].z * pos2[i].y; C.zz() += pos1[i].z * pos2[i].z; } // build the "overlap" matrix, whose eigenvectors are stationary // points of the RMSD in the space of rotations S[0][0] = C.xx() + C.yy() + C.zz(); S[1][0] = C.yz() - C.zy(); S[0][1] = S[1][0]; S[2][0] = - C.xz() + C.zx() ; S[0][2] = S[2][0]; S[3][0] = C.xy() - C.yx(); S[0][3] = S[3][0]; S[1][1] = C.xx() - C.yy() - C.zz(); S[2][1] = C.xy() + C.yx(); S[1][2] = S[2][1]; S[3][1] = C.xz() + C.zx(); S[1][3] = S[3][1]; S[2][2] = - C.xx() + C.yy() - C.zz(); S[3][2] = C.yz() + C.zy(); S[2][3] = S[3][2]; S[3][3] = - C.xx() - C.yy() + C.zz(); // if (cvm::debug()) { // for (size_t i = 0; i < 4; i++) { // std::string line (""); // for (size_t j = 0; j < 4; j++) { // line += std::string (" S["+cvm::to_str (i)+ // "]["+cvm::to_str (j)+"] ="+cvm::to_str (S[i][j])); // } // cvm::log (line+"\n"); // } // } } void colvarmodule::rotation::diagonalize_matrix (matrix2d &S, cvm::real S_eigval[4], matrix2d &S_eigvec) { // diagonalize int jac_nrot = 0; jacobi (S, 4, S_eigval, S_eigvec, &jac_nrot); eigsrt (S_eigval, S_eigvec, 4); // jacobi saves eigenvectors by columns transpose (S_eigvec, 4); // normalize eigenvectors for (size_t ie = 0; ie < 4; ie++) { cvm::real norm2 = 0.0; for (size_t i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2)); cvm::real const norm = std::sqrt (norm2); for (size_t i = 0; i < 4; i++) S_eigvec[ie][i] /= norm; } } // Calculate the rotation, plus its derivatives void colvarmodule::rotation::calc_optimal_rotation (std::vector const &pos1, std::vector const &pos2) { matrix2d S; matrix2d S_backup; cvm::real S_eigval[4]; matrix2d S_eigvec; // if (cvm::debug()) { // cvm::atom_pos cog1 (0.0, 0.0, 0.0); // for (size_t i = 0; i < pos1.size(); i++) { // cog1 += pos1[i]; // } // cog1 /= cvm::real (pos1.size()); // cvm::atom_pos cog2 (0.0, 0.0, 0.0); // for (size_t i = 0; i < pos2.size(); i++) { // cog2 += pos2[i]; // } // cog2 /= cvm::real (pos1.size()); // cvm::log ("calc_optimal_rotation: centers of geometry are: "+ // cvm::to_str (cog1, cvm::cv_width, cvm::cv_prec)+ // " and "+cvm::to_str (cog2, cvm::cv_width, cvm::cv_prec)+".\n"); // } build_matrix (pos1, pos2, S); S_backup = S; if (cvm::debug()) { if (b_debug_gradients) { cvm::log ("S = "+cvm::to_str (cvm::to_str (S_backup), cvm::cv_width, cvm::cv_prec)+"\n"); } } diagonalize_matrix (S, S_eigval, S_eigvec); // eigenvalues and eigenvectors cvm::real const &L0 = S_eigval[0]; cvm::real const &L1 = S_eigval[1]; cvm::real const &L2 = S_eigval[2]; cvm::real const &L3 = S_eigval[3]; cvm::real const *Q0 = S_eigvec[0]; cvm::real const *Q1 = S_eigvec[1]; cvm::real const *Q2 = S_eigvec[2]; cvm::real const *Q3 = S_eigvec[3]; lambda = L0; q = cvm::quaternion (Q0); if (q_old.norm2() > 0.0) { q.match (q_old); if (q_old.inner (q) < (1.0 - crossing_threshold)) { cvm::log ("Warning: discontinuous rotation!\n"); } } q_old = q; if (cvm::debug()) { if (b_debug_gradients) { cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+ ", Q0 = "+cvm::to_str (cvm::quaternion (Q0), cvm::cv_width, cvm::cv_prec)+ ", Q0*Q0 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q0)), cvm::cv_width, cvm::cv_prec)+ "\n"); cvm::log ("L1 = "+cvm::to_str (L1, cvm::cv_width, cvm::cv_prec)+ ", Q1 = "+cvm::to_str (cvm::quaternion (Q1), cvm::cv_width, cvm::cv_prec)+ ", Q0*Q1 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q1)), cvm::cv_width, cvm::cv_prec)+ "\n"); cvm::log ("L2 = "+cvm::to_str (L2, cvm::cv_width, cvm::cv_prec)+ ", Q2 = "+cvm::to_str (cvm::quaternion (Q2), cvm::cv_width, cvm::cv_prec)+ ", Q0*Q2 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q2)), cvm::cv_width, cvm::cv_prec)+ "\n"); cvm::log ("L3 = "+cvm::to_str (L3, cvm::cv_width, cvm::cv_prec)+ ", Q3 = "+cvm::to_str (cvm::quaternion (Q3), cvm::cv_width, cvm::cv_prec)+ ", Q0*Q3 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q3)), cvm::cv_width, cvm::cv_prec)+ "\n"); } } // calculate derivatives of L0 and Q0 with respect to each atom in // either group; note: if dS_1 is a null vector, nothing will be // calculated for (size_t ia = 0; ia < dS_1.size(); ia++) { cvm::real const &a2x = pos2[ia].x; cvm::real const &a2y = pos2[ia].y; cvm::real const &a2z = pos2[ia].z; matrix2d &ds_1 = dS_1[ia]; // derivative of the S matrix ds_1.reset(); ds_1[0][0] = cvm::rvector ( a2x, a2y, a2z); ds_1[1][0] = cvm::rvector ( 0.0, a2z, -a2y); ds_1[0][1] = ds_1[1][0]; ds_1[2][0] = cvm::rvector (-a2z, 0.0, a2x); ds_1[0][2] = ds_1[2][0]; ds_1[3][0] = cvm::rvector ( a2y, -a2x, 0.0); ds_1[0][3] = ds_1[3][0]; ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z); ds_1[2][1] = cvm::rvector ( a2y, a2x, 0.0); ds_1[1][2] = ds_1[2][1]; ds_1[3][1] = cvm::rvector ( a2z, 0.0, a2x); ds_1[1][3] = ds_1[3][1]; ds_1[2][2] = cvm::rvector (-a2x, a2y, -a2z); ds_1[3][2] = cvm::rvector ( 0.0, a2z, a2y); ds_1[2][3] = ds_1[3][2]; ds_1[3][3] = cvm::rvector (-a2x, -a2y, a2z); cvm::rvector &dl0_1 = dL0_1[ia]; vector1d &dq0_1 = dQ0_1[ia]; // matrix multiplications; derivatives of L_0 and Q_0 are // calculated using Hellmann-Feynman theorem (i.e. exploiting the // fact that the eigenvectors Q_i form an orthonormal basis) dl0_1.reset(); for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dl0_1 += Q0[i] * ds_1[i][j] * Q0[j]; } } dq0_1.reset(); for (size_t p = 0; p < 4; p++) { for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dq0_1[p] += (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] + (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] + (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p]; } } } } // do the same for the second group for (size_t ia = 0; ia < dS_2.size(); ia++) { cvm::real const &a1x = pos1[ia].x; cvm::real const &a1y = pos1[ia].y; cvm::real const &a1z = pos1[ia].z; matrix2d &ds_2 = dS_2[ia]; ds_2.reset(); ds_2[0][0] = cvm::rvector ( a1x, a1y, a1z); ds_2[1][0] = cvm::rvector ( 0.0, -a1z, a1y); ds_2[0][1] = ds_2[1][0]; ds_2[2][0] = cvm::rvector ( a1z, 0.0, -a1x); ds_2[0][2] = ds_2[2][0]; ds_2[3][0] = cvm::rvector (-a1y, a1x, 0.0); ds_2[0][3] = ds_2[3][0]; ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z); ds_2[2][1] = cvm::rvector ( a1y, a1x, 0.0); ds_2[1][2] = ds_2[2][1]; ds_2[3][1] = cvm::rvector ( a1z, 0.0, a1x); ds_2[1][3] = ds_2[3][1]; ds_2[2][2] = cvm::rvector (-a1x, a1y, -a1z); ds_2[3][2] = cvm::rvector ( 0.0, a1z, a1y); ds_2[2][3] = ds_2[3][2]; ds_2[3][3] = cvm::rvector (-a1x, -a1y, a1z); cvm::rvector &dl0_2 = dL0_2[ia]; vector1d &dq0_2 = dQ0_2[ia]; dl0_2.reset(); for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dl0_2 += Q0[i] * ds_2[i][j] * Q0[j]; } } dq0_2.reset(); for (size_t p = 0; p < 4; p++) { for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dq0_2[p] += (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] + (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] + (Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p]; } } } if (cvm::debug()) { if (b_debug_gradients) { matrix2d S_new; cvm::real S_new_eigval[4]; matrix2d S_new_eigvec; // make an infitesimal move along each cartesian coordinate of // this atom, and solve again the eigenvector problem for (size_t comp = 0; comp < 3; comp++) { S_new = S_backup; // diagonalize the new overlap matrix for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { S_new[i][j] += colvarmodule::debug_gradients_step_size * ds_2[i][j][comp]; } } // cvm::log ("S_new = "+cvm::to_str (cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n"); diagonalize_matrix (S_new, S_new_eigval, S_new_eigvec); cvm::real const &L0_new = S_new_eigval[0]; cvm::real const *Q0_new = S_new_eigvec[0]; cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size; cvm::quaternion const q0 (Q0); cvm::quaternion const DQ0 (dq0_2[0][comp] * colvarmodule::debug_gradients_step_size, dq0_2[1][comp] * colvarmodule::debug_gradients_step_size, dq0_2[2][comp] * colvarmodule::debug_gradients_step_size, dq0_2[3][comp] * colvarmodule::debug_gradients_step_size); cvm::log ( "|(l_0+dl_0) - l_0^new|/l_0 = "+ cvm::to_str (std::fabs (L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+ ", |(q_0+dq_0) - q_0^new| = "+ cvm::to_str ((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+ "\n"); } } } } } // Numerical Recipes routine for diagonalization #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau); \ a[k][l]=h+s*(g-h*tau); void jacobi(cvm::real **a, int n, cvm::real d[], cvm::real **v, int *nrot) { int j,iq,ip,i; cvm::real tresh,theta,tau,t,sm,s,h,g,c; std::vector b (n, 0.0); std::vector z (n, 0.0); for (ip=0;ip 4 && (cvm::real)(std::fabs(d[ip])+g) == (cvm::real)std::fabs(d[ip]) && (cvm::real)(std::fabs(d[iq])+g) == (cvm::real)std::fabs(d[iq])) a[ip][iq]=0.0; else if (std::fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h)) t=(a[ip][iq])/h; else { theta=0.5*h/(a[ip][iq]); t=1.0/(std::fabs(theta)+std::sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } c=1.0/std::sqrt(1+t*t); s=t*c; tau=s/(1.0+c); h=t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq]=0.0; for (j=0;j<=ip-1;j++) { ROTATE(a,j,ip,j,iq) } for (j=ip+1;j<=iq-1;j++) { ROTATE(a,ip,j,j,iq) } for (j=iq+1;j= p) p=d[k=j]; if (k != i) { d[k]=d[i]; d[i]=p; for (j=0;j #include #include #include #include #include #include #include class colvarparse; class colvar; class colvarbias; class colvarproxy; /// \brief Collective variables module (main class) /// /// Class to control the collective variables calculation. An object /// (usually one) of this class is spawned from the MD program, /// containing all i/o routines and general interface. /// /// At initialization, the colvarmodule object creates a proxy object /// to provide a transparent interface between the MD program and the /// child objects class colvarmodule { private: /// Impossible to initialize the main object without arguments colvarmodule(); public: friend class colvarproxy; /// Defining an abstract real number allows to switch precision typedef double real; /// Residue identifier typedef int residue_id; class rvector; template class vector1d; template class matrix2d; class quaternion; class rotation; /// \brief Atom position (different type name from rvector, to make /// possible future PBC-transparent implementations) typedef rvector atom_pos; /// \brief 3x3 matrix of real numbers class rmatrix; // allow these classes to access protected data class atom; class atom_group; friend class atom; friend class atom_group; typedef std::vector::iterator atom_iter; typedef std::vector::const_iterator atom_const_iter; /// Current step number static size_t it; /// Starting step number for this run static size_t it_restart; /// Return the current step number from the beginning of this run static inline size_t step_relative() { return it - it_restart; } /// Return the current step number from the beginning of the whole /// calculation static inline size_t step_absolute() { return it; } /// If true, get it_restart from the state file; if set to false, /// the MD program is providing it bool it_restart_from_state_file; /// \brief Finite difference step size (if there is no dynamics, or /// if gradients need to be tested independently from the size of /// dt) static real debug_gradients_step_size; /// Prefix for all output files for this run static std::string output_prefix; /// Prefix for files from a previous run (including restart/output) static std::string input_prefix; /// input restart file name (determined from input_prefix) static std::string restart_in_name; /// Array of collective variables static std::vector colvars; /* TODO: implement named CVCs /// Array of named (reusable) collective variable components static std::vector cvcs; /// Named cvcs register themselves at initialization time inline void register_cvc (cvc *p) { cvcs.push_back(p); } */ /// Array of collective variable biases static std::vector biases; /// \brief Number of ABF biases initialized (in normal conditions /// should be 1) static size_t n_abf_biases; /// \brief Number of metadynamics biases initialized (in normal /// conditions should be 1) static size_t n_meta_biases; /// \brief Number of harmonic biases initialized (no limit on the /// number) static size_t n_harm_biases; /// \brief Number of histograms initialized (no limit on the /// number) static size_t n_histo_biases; /// \brief Whether debug output should be enabled (compile-time option) static inline bool debug() { return COLVARS_DEBUG; } /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name colvarmodule (char const *config_name, colvarproxy *proxy_in); /// Destructor ~colvarmodule(); /// Initialize collective variables void init_colvars (std::string const &conf); /// Initialize collective variable biases void init_biases (std::string const &conf); /// Load new configuration - force constant and/or centers only void change_configuration(std::string const &name, std::string const &conf); /// Calculate change in energy from using alternate configuration real energy_difference(std::string const &name, std::string const &conf); /// Calculate collective variables and biases void calc(); /// Read the input restart file std::istream & read_restart (std::istream &is); /// Write the output restart file std::ostream & write_restart (std::ostream &os); /// Write all output files (called by the proxy) void write_output_files(); /// \brief Call colvarproxy::backup_file() static void backup_file (char const *filename); /// Perform analysis void analyze(); /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) bool read_traj (char const *traj_filename, size_t traj_read_begin, size_t traj_read_end); /// Get the pointer of a colvar from its name (returns NULL if not found) static colvar * colvar_p (std::string const &name); /// Quick conversion of an object to a string template static std::string to_str (T const &x, size_t const &width = 0, size_t const &prec = 0); /// Quick conversion of a vector of objects to a string template static std::string to_str (std::vector const &x, size_t const &width = 0, size_t const &prec = 0); /// Reduce the number of characters in a string static inline std::string wrap_string (std::string const &s, size_t const &nchars) { if (!s.size()) return std::string (nchars, ' '); else return ( (s.size() <= size_t (nchars)) ? (s+std::string (nchars-s.size(), ' ')) : (std::string (s, 0, nchars)) ); } /// Number of characters to represent a time step static size_t const it_width; /// Number of digits to represent a collective variables value(s) static size_t const cv_prec; /// Number of characters to represent a collective variables value(s) static size_t const cv_width; /// Number of digits to represent the collective variables energy static size_t const en_prec; /// Number of characters to represent the collective variables energy static size_t const en_width; /// Line separator in the log output static std::string const line_marker; // proxy functions /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) static real unit_angstrom(); /// \brief Boltmann constant static real boltzmann(); /// \brief Temperature of the simulation (K) static real temperature(); /// \brief Time step of MD integrator (fs) static real dt(); /// Request calculation of system force from MD engine static void request_system_force(); /// Print a message to the main log static void log (std::string const &message); /// Print a message to the main log and exit with error code static void fatal_error (std::string const &message); /// Print a message to the main log and exit normally static void exit (std::string const &message); /// \brief Get the distance between two atomic positions with pbcs handled /// correctly static rvector position_distance (atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the square distance between two positions (with /// periodic boundary conditions handled transparently) /// /// Note: in the case of periodic boundary conditions, this provides /// an analytical square distance (while taking the square of /// position_distance() would produce leads to a cusp) static real position_dist2 (atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image /// \param ref_pos (optional) The reference position static void select_closest_image (atom_pos &pos, atom_pos const &ref_pos); /// \brief Perform select_closest_image() on a set of atomic positions /// /// After that, distance vectors can then be calculated directly, /// without using position_distance() static void select_closest_images (std::vector &pos, atom_pos const &ref_pos); /// \brief Create atoms from a file \param filename name of the file /// (usually a PDB) \param atoms array of the atoms to be allocated /// \param pdb_field (optiona) if "filename" is a PDB file, use this /// field to determine which are the atoms to be set static void load_atoms (char const *filename, std::vector &atoms, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from a file /// (usually a PDB); the number of atoms in "filename" must match /// the number of elements in "pos" static void load_coords (char const *filename, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value = 0.0); /// Frequency for collective variables trajectory output static size_t cv_traj_freq; /// \brief True if only analysis is performed and not a run static bool b_analysis; /// Frequency for saving output restarts static size_t restart_out_freq; /// Output restart file name std::string restart_out_name; /// Pseudo-random number with Gaussian distribution static real rand_gaussian (void); protected: /// Configuration file std::ifstream config_s; /// Configuration file parser object colvarparse *parse; /// Name of the trajectory file std::string cv_traj_name; /// Collective variables output trajectory file std::ofstream cv_traj_os; /// Appending to the existing trajectory file? bool cv_traj_append; /// Output restart file std::ofstream restart_out_os; /// \brief Pointer to the proxy object, used to retrieve atomic data /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; /// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu static size_t depth; public: /// Increase the depth (number of indentations in the output) static void increase_depth(); /// Decrease the depth (number of indentations in the output) static void decrease_depth(); }; /// Shorthand for the frequently used type prefix typedef colvarmodule cvm; #include "colvartypes.h" std::ostream & operator << (std::ostream &os, cvm::rvector const &v); std::istream & operator >> (std::istream &is, cvm::rvector &v); template std::string cvm::to_str (T const &x, size_t const &width, size_t const &prec) { std::ostringstream os; if (width) os.width (width); if (prec) { os.setf (std::ios::scientific, std::ios::floatfield); os.precision (prec); } os << x; return os.str(); } template std::string cvm::to_str (std::vector const &x, size_t const &width, size_t const &prec) { if (!x.size()) return std::string (""); std::ostringstream os; if (prec) { os.setf (std::ios::scientific, std::ios::floatfield); } os << "{ "; if (width) os.width (width); if (prec) os.precision (prec); os << x[0]; for (size_t i = 1; i < x.size(); i++) { os << ", "; if (width) os.width (width); if (prec) os.precision (prec); os << x[i]; } os << " }"; return os.str(); } #include "colvarproxy.h" inline cvm::real cvm::unit_angstrom() { return proxy->unit_angstrom(); } inline cvm::real cvm::boltzmann() { return proxy->boltzmann(); } inline cvm::real cvm::temperature() { return proxy->temperature(); } inline cvm::real cvm::dt() { return proxy->dt(); } inline void cvm::request_system_force() { proxy->request_system_force (true); } inline void cvm::select_closest_image (atom_pos &pos, atom_pos const &ref_pos) { proxy->select_closest_image (pos, ref_pos); } inline void cvm::select_closest_images (std::vector &pos, atom_pos const &ref_pos) { proxy->select_closest_images (pos, ref_pos); } inline cvm::rvector cvm::position_distance (atom_pos const &pos1, atom_pos const &pos2) { return proxy->position_distance (pos1, pos2); } inline cvm::real cvm::position_dist2 (cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { return proxy->position_dist2 (pos1, pos2); } inline void cvm::load_atoms (char const *file_name, std::vector &atoms, std::string const &pdb_field, double const pdb_field_value) { proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value); } inline void cvm::load_coords (char const *file_name, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value) { proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value); } inline void cvm::backup_file (char const *filename) { proxy->backup_file (filename); } inline cvm::real cvm::rand_gaussian (void) { return proxy->rand_gaussian(); } #endif // Emacs // Local Variables: // mode: C++ // End: