diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 07d8254475..6194424212 100644 Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ diff --git a/doc/src/fix_colvars.txt b/doc/src/fix_colvars.txt index d9eb9ac5aa..304273d451 100644 --- a/doc/src/fix_colvars.txt +++ b/doc/src/fix_colvars.txt @@ -98,6 +98,16 @@ fix to add the energy change from the biasing force added by the fix to the system's potential energy as part of "thermodynamic output"_thermo_style.html. +The {fix_modify configfile } option allows to add settings +from an additional config file to the colvars module. This option can +only be used, after the system has been initialized with a "run"_run.html +command. + +The {fix_modify config } option allows to add settings +from inline strings. Those have to fit on a single line when enclosed +in a pair of double quotes ("), or can span multiple lines when bracketed +by a pair of triple double quotes (""", like python embedded documentation). + This fix computes a global scalar which can be accessed by various "output commands"_Howto_output.html. The scalar is the cumulative energy change due to this fix. The scalar value calculated by this diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 8b28eaa0df..723a54ad39 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -24,12 +24,13 @@ colvar::colvar() { runave_os = NULL; - prev_timestep = -1; + prev_timestep = -1L; after_restart = false; kinetic_energy = 0.0; potential_energy = 0.0; - init_cv_requires(); + description = "uninitialized colvar"; + init_dependencies(); } @@ -193,7 +194,7 @@ int colvar::init(std::string const &conf) { bool homogeneous = is_enabled(f_cv_linear); for (i = 0; i < cvcs.size(); i++) { - if ((std::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { + if ((cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { homogeneous = false; } } @@ -224,7 +225,7 @@ int colvar::init(std::string const &conf) // Allow scripted/custom functions to be defined as periodic if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) { if (get_keyval(conf, "period", period, 0.)) { - set_enabled(f_cv_periodic, true); + enable(f_cv_periodic); get_keyval(conf, "wrapAround", wrap_center, 0.); } } @@ -471,7 +472,7 @@ int colvar::init_grid_parameters(std::string const &conf) if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) { cvm::log("Reading legacy options lowerWall and lowerWallConstant: " - "consider using a harmonicWalls restraint.\n"); + "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n"); lower_wall.type(value()); if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) { cvm::log("Warning: lowerWall will need to be " @@ -485,7 +486,7 @@ int colvar::init_grid_parameters(std::string const &conf) if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) { cvm::log("Reading legacy options upperWall and upperWallConstant: " - "consider using a harmonicWalls restraint.\n"); + "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n"); upper_wall.type(value()); if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) { cvm::log("Warning: upperWall will need to be " @@ -562,13 +563,13 @@ int colvar::init_extended_Lagrangian(std::string const &conf) get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false); if (is_enabled(f_cv_extended_Lagrangian)) { - cvm::real temp, tolerance, period; + cvm::real temp, tolerance, extended_period; cvm::log("Enabling the extended Lagrangian term for colvar \""+ this->name+"\".\n"); - xr.type(value()); - vr.type(value()); + x_ext.type(value()); + v_ext.type(value()); fr.type(value()); const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); @@ -590,11 +591,11 @@ int colvar::init_extended_Lagrangian(std::string const &conf) ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2"); - get_keyval(conf, "extendedTimeConstant", period, 200.0); - if (period <= 0.0) { + get_keyval(conf, "extendedTimeConstant", extended_period, 200.0); + if (extended_period <= 0.0) { cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR); } - ext_mass = (cvm::boltzmann() * temp * period * period) + ext_mass = (cvm::boltzmann() * temp * extended_period * extended_period) / (4.0 * PI * PI * tolerance * tolerance); cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)"); @@ -615,7 +616,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf) enable(f_cv_Langevin); ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs // Adjust Langevin sigma for slow time step if time_step_factor != 1 - ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor))); + ext_sigma = cvm::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor))); } } @@ -761,6 +762,8 @@ int colvar::init_components(std::string const &conf) "weighted by inverse power", "distanceInv"); error_code |= init_components_type(conf, "N1xN2-long vector " "of pairwise distances", "distancePairs"); + error_code |= init_components_type(conf, "dipole magnitude", + "dipoleMagnitude"); error_code |= init_components_type(conf, "coordination " "number", "coordNum"); error_code |= init_components_type(conf, "self-coordination " @@ -831,22 +834,25 @@ void colvar::build_atom_list(void) for (size_t i = 0; i < cvcs.size(); i++) { for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { - cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + cvm::atom_group const &ag = *(cvcs[i]->atom_groups[j]); for (size_t k = 0; k < ag.size(); k++) { temp_id_list.push_back(ag[k].id); } + if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) { + cvm::atom_group const &fg = *(ag.fitting_group); + for (size_t k = 0; k < fg.size(); k++) { + temp_id_list.push_back(fg[k].id); + } + } } } temp_id_list.sort(); temp_id_list.unique(); - // atom_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); - unsigned int id_i = 0; std::list::iterator li; for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) { - atom_ids[id_i] = *li; - id_i++; + atom_ids.push_back(*li); } temp_id_list.clear(); @@ -934,16 +940,153 @@ int colvar::parse_analysis(std::string const &conf) } -void colvar::setup() { - // loop over all components to reset masses of all groups - 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.setup(); - atoms.reset_mass(name,i,ig); - atoms.read_positions(); +int colvar::init_dependencies() { + size_t i; + if (features().size() == 0) { + for (i = 0; i < f_cv_ntot; i++) { + modify_features().push_back(new feature); + } + + init_feature(f_cv_active, "active", f_type_dynamic); + // Do not require f_cvc_active in children, as some components may be disabled + // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom + require_feature_alt(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function); + + init_feature(f_cv_awake, "awake", f_type_static); + require_feature_self(f_cv_awake, f_cv_active); + + init_feature(f_cv_gradient, "gradient", f_type_dynamic); + require_feature_children(f_cv_gradient, f_cvc_gradient); + + init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic); + require_feature_self(f_cv_collect_gradient, f_cv_gradient); + require_feature_self(f_cv_collect_gradient, f_cv_scalar); + // The following exlusion could be lifted by implementing the feature + exclude_feature_self(f_cv_collect_gradient, f_cv_scripted); + require_feature_children(f_cv_collect_gradient, f_cvc_explicit_gradient); + + init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic); + + // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly + init_feature(f_cv_total_force, "total force", f_type_dynamic); + require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc); + + // Deps for explicit total force calculation + init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic); + require_feature_self(f_cv_total_force_calc, f_cv_scalar); + require_feature_self(f_cv_total_force_calc, f_cv_linear); + require_feature_children(f_cv_total_force_calc, f_cvc_inv_gradient); + require_feature_self(f_cv_total_force_calc, f_cv_Jacobian); + + init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic); + require_feature_self(f_cv_Jacobian, f_cv_scalar); + require_feature_self(f_cv_Jacobian, f_cv_linear); + require_feature_children(f_cv_Jacobian, f_cvc_Jacobian); + + init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user); + require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated + + init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user); + require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar); + require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient); + + init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user); + require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian); + + init_feature(f_cv_linear, "linear", f_type_static); + + init_feature(f_cv_scalar, "scalar", f_type_static); + + init_feature(f_cv_output_energy, "output energy", f_type_user); + + init_feature(f_cv_output_value, "output value", f_type_user); + + init_feature(f_cv_output_velocity, "output velocity", f_type_user); + require_feature_self(f_cv_output_velocity, f_cv_fdiff_velocity); + + init_feature(f_cv_output_applied_force, "output applied force", f_type_user); + + init_feature(f_cv_output_total_force, "output total force", f_type_user); + require_feature_self(f_cv_output_total_force, f_cv_total_force); + + init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user); + require_feature_self(f_cv_subtract_applied_force, f_cv_total_force); + + init_feature(f_cv_lower_boundary, "lower boundary", f_type_user); + require_feature_self(f_cv_lower_boundary, f_cv_scalar); + + init_feature(f_cv_upper_boundary, "upper boundary", f_type_user); + require_feature_self(f_cv_upper_boundary, f_cv_scalar); + + init_feature(f_cv_grid, "grid", f_type_dynamic); + require_feature_self(f_cv_grid, f_cv_lower_boundary); + require_feature_self(f_cv_grid, f_cv_upper_boundary); + + init_feature(f_cv_runave, "running average", f_type_user); + + init_feature(f_cv_corrfunc, "correlation function", f_type_user); + + init_feature(f_cv_scripted, "scripted", f_type_user); + + init_feature(f_cv_custom_function, "custom function", f_type_user); + exclude_feature_self(f_cv_custom_function, f_cv_scripted); + + init_feature(f_cv_periodic, "periodic", f_type_static); + require_feature_self(f_cv_periodic, f_cv_scalar); + init_feature(f_cv_scalar, "scalar", f_type_static); + init_feature(f_cv_linear, "linear", f_type_static); + init_feature(f_cv_homogeneous, "homogeneous", f_type_static); + + // because total forces are obtained from the previous time step, + // we cannot (currently) have colvar values and total forces for the same timestep + init_feature(f_cv_multiple_ts, "multiple timestep colvar", f_type_static); + exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc); + + // check that everything is initialized + for (i = 0; i < colvardeps::f_cv_ntot; i++) { + if (is_not_set(i)) { + cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description); + } } } + + // Initialize feature_states for each instance + feature_states.reserve(f_cv_ntot); + for (i = 0; i < f_cv_ntot; i++) { + feature_states.push_back(feature_state(true, false)); + // Most features are available, so we set them so + // and list exceptions below + } + + feature_states[f_cv_fdiff_velocity].available = + cvm::main()->proxy->simulation_running(); + + return COLVARS_OK; +} + + +void colvar::setup() +{ + // loop over all components to update masses and charges of all groups + 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->setup(); + atoms->print_properties(name, i, ig); + atoms->read_positions(); + } + } +} + + +std::vector > colvar::get_atom_lists() +{ + std::vector > lists; + for (size_t i = 0; i < cvcs.size(); i++) { + std::vector > li = cvcs[i]->get_atom_lists(); + lists.insert(lists.end(), li.begin(), li.end()); + } + return lists; } @@ -953,8 +1096,8 @@ colvar::~colvar() // because the children are cvcs and will be deleted // just below -// Clear references to this colvar's cvcs as children -// for dependency purposes + // Clear references to this colvar's cvcs as children + // for dependency purposes remove_all_children(); for (std::vector::reverse_iterator ci = cvcs.rbegin(); @@ -1231,7 +1374,6 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs) int colvar::collect_cvc_gradients() { size_t i; - if (is_enabled(f_cv_collect_gradient)) { // Collect the atomic gradients inside colvar object for (unsigned int a = 0; a < atomic_gradients.size(); a++) { @@ -1239,34 +1381,7 @@ int colvar::collect_cvc_gradients() } for (i = 0; i < cvcs.size(); i++) { if (!cvcs[i]->is_enabled()) continue; - // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx - cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) * - cvm::integer_power((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); - - for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { - - cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); - - // If necessary, apply inverse rotation to get atomic - // gradient in the laboratory frame - if (ag.b_rotate) { - cvm::rotation const rot_inv = ag.rot.inverse(); - - for (size_t k = 0; k < ag.size(); k++) { - size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - ag[k].id) - atom_ids.begin(); - atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad); - } - - } else { - - for (size_t k = 0; k < ag.size(); k++) { - size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - ag[k].id) - atom_ids.begin(); - atomic_gradients[a] += coeff * ag[k].grad; - } - } - } + cvcs[i]->collect_gradients(atom_ids, atomic_gradients); } } return COLVARS_OK; @@ -1391,20 +1506,20 @@ int colvar::calc_colvar_properties() // initialize the restraint center in the first step to the value // just calculated from the cvcs if (cvm::step_relative() == 0 && !after_restart) { - xr = x; - vr.reset(); // (already 0; added for clarity) + x_ext = x; + v_ext.reset(); // (already 0; added for clarity) } // Special case of a repeated timestep (eg. multiple NAMD "run" statements) // revert values of the extended coordinate and velocity prior to latest integration - if (cvm::step_relative() == prev_timestep) { - xr = prev_xr; - vr = prev_vr; + if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) { + x_ext = prev_x_ext; + v_ext = prev_v_ext; } // report the restraint center as "value" - x_reported = xr; - v_reported = vr; + x_reported = x_ext; + v_reported = v_ext; // the "total force" with the extended Lagrangian is // calculated in update_forces_energy() below @@ -1458,77 +1573,86 @@ cvm::real colvar::update_forces_energy() // extended variable if there is one if (is_enabled(f_cv_extended_Lagrangian)) { - if (cvm::debug()) { - cvm::log("Updating extended-Lagrangian degree of freedom.\n"); - } - - if (prev_timestep > -1) { - // Keep track of slow timestep to integrate MTS colvars - // the colvar checks the interval after waking up twice - int n_timesteps = cvm::step_relative() - prev_timestep; - if (n_timesteps != 0 && n_timesteps != time_step_factor) { - cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " + - cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) + - " steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " + - cvm::to_str(cvm::step_relative()) + ").\n" + - "Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n"); - return 0.; + if (cvm::proxy->simulation_running()) { + // Only integrate the extended equations of motion in running MD simulations + if (cvm::debug()) { + cvm::log("Updating extended-Lagrangian degree of freedom.\n"); } - } - // Integrate with slow timestep (if time_step_factor != 1) - cvm::real dt = cvm::dt() * cvm::real(time_step_factor); + if (prev_timestep > -1) { + // Keep track of slow timestep to integrate MTS colvars + // the colvar checks the interval after waking up twice + int n_timesteps = cvm::step_relative() - prev_timestep; + if (n_timesteps != 0 && n_timesteps != time_step_factor) { + cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " + + cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) + + " steps at timestep " + cvm::to_str(cvm::step_absolute()) + " (relative step: " + + cvm::to_str(cvm::step_relative()) + ").\n" + + "Make sure that this colvar is requested by biases at multiples of timeStepFactor.\n"); + return 0.; + } + } - colvarvalue f_ext(fr.type()); // force acting on the extended variable - f_ext.reset(); + // Integrate with slow timestep (if time_step_factor != 1) + cvm::real dt = cvm::dt() * cvm::real(time_step_factor); - // the total force is applied to the fictitious mass, while the - // atoms only feel the harmonic force + wall force - // fr: bias force on extended variable (without harmonic spring), for output in trajectory - // f_ext: total force on extended variable (including harmonic spring) - // f: - initially, external biasing force - // - after this code block, colvar force to be applied to atomic coordinates - // ie. spring force (fb_actual will be added just below) - fr = f; - // External force has been scaled for a 1-timestep impulse, scale it back because we will - // integrate it with the colvar's own timestep factor - f_ext = f / cvm::real(time_step_factor); - f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); - f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x); - // Coupling force is a slow force, to be applied to atomic coords impulse-style - f *= cvm::real(time_step_factor); + colvarvalue f_ext(fr.type()); // force acting on the extended variable + f_ext.reset(); - if (is_enabled(f_cv_subtract_applied_force)) { - // Report a "system" force without the biases on this colvar - // that is, just the spring force - ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); + // the total force is applied to the fictitious mass, while the + // atoms only feel the harmonic force + wall force + // fr: bias force on extended variable (without harmonic spring), for output in trajectory + // f_ext: total force on extended variable (including harmonic spring) + // f: - initially, external biasing force + // - after this code block, colvar force to be applied to atomic coordinates + // ie. spring force (fb_actual will be added just below) + fr = f; + // External force has been scaled for a 1-timestep impulse, scale it back because we will + // integrate it with the colvar's own timestep factor + f_ext = f / cvm::real(time_step_factor); + f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x); + f = (-0.5 * ext_force_k) * this->dist2_rgrad(x_ext, x); + // Coupling force is a slow force, to be applied to atomic coords impulse-style + f *= cvm::real(time_step_factor); + + if (is_enabled(f_cv_subtract_applied_force)) { + // Report a "system" force without the biases on this colvar + // that is, just the spring force + ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x); + } else { + // The total force acting on the extended variable is f_ext + // This will be used in the next timestep + ft_reported = f_ext; + } + + // backup in case we need to revert this integration timestep + // if the same MD timestep is re-run + prev_x_ext = x_ext; + prev_v_ext = v_ext; + + // leapfrog: starting from x_i, f_i, v_(i-1/2) + v_ext += (0.5 * dt) * f_ext / ext_mass; + // Because of leapfrog, kinetic energy at time i is approximate + kinetic_energy = 0.5 * ext_mass * v_ext * v_ext; + potential_energy = 0.5 * ext_force_k * this->dist2(x_ext, x); + // leap to v_(i+1/2) + if (is_enabled(f_cv_Langevin)) { + v_ext -= dt * ext_gamma * v_ext; + colvarvalue rnd(x); + rnd.set_random(); + v_ext += dt * ext_sigma * rnd / ext_mass; + } + v_ext += (0.5 * dt) * f_ext / ext_mass; + x_ext += dt * v_ext; + x_ext.apply_constraints(); + this->wrap(x_ext); } else { - // The total force acting on the extended variable is f_ext - // This will be used in the next timestep - ft_reported = f_ext; + // If this is a postprocessing run (eg. in VMD), the extended DOF + // is equal to the actual coordinate + x_ext = x; } - - // backup in case we need to revert this integration timestep - // if the same MD timestep is re-run - prev_xr = xr; - prev_vr = vr; - - // leapfrog: starting from x_i, f_i, v_(i-1/2) - vr += (0.5 * dt) * f_ext / 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 (is_enabled(f_cv_Langevin)) { - vr -= dt * ext_gamma * vr; - colvarvalue rnd(x); - rnd.set_random(); - vr += dt * ext_sigma * rnd / ext_mass; - } - vr += (0.5 * dt) * f_ext / ext_mass; - xr += dt * vr; - xr.apply_constraints(); - this->wrap(xr); + // Report extended value + x_reported = x_ext; } // Now adding the force on the actual colvar (for those biases that @@ -1730,7 +1854,7 @@ int colvar::update_cvc_config(std::vector const &confs) // ******************** METRIC FUNCTIONS ******************** -// Use the metrics defined by \link cvc \endlink objects +// Use the metrics defined by \link colvar::cvc \endlink objects bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const @@ -1742,7 +1866,7 @@ bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) c } if (period > 0.0) { - if ( ((std::sqrt(this->dist2(lb, ub))) / this->width) + if ( ((cvm::sqrt(this->dist2(lb, ub))) / this->width) < 1.0E-10 ) { return true; } @@ -1792,21 +1916,21 @@ colvarvalue colvar::dist2_rgrad(colvarvalue const &x1, } } -void colvar::wrap(colvarvalue &x) const + +void colvar::wrap(colvarvalue &x_unwrapped) const { - if ( !is_enabled(f_cv_periodic) ) { + if (!is_enabled(f_cv_periodic)) { return; } if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) { // Scripted functions do their own wrapping, as cvcs might not be periodic - cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5); - x.real_value -= shift * period; + cvm::real shift = cvm::floor((x_unwrapped.real_value - wrap_center) / + period + 0.5); + x_unwrapped.real_value -= shift * period; } else { - cvcs[0]->wrap(x); + cvcs[0]->wrap(x_unwrapped); } - - return; } @@ -1852,15 +1976,15 @@ std::istream & colvar::read_restart(std::istream &is) if (is_enabled(f_cv_extended_Lagrangian)) { - if ( !(get_keyval(conf, "extended_x", xr, + if ( !(get_keyval(conf, "extended_x", x_ext, colvarvalue(x.type()), colvarparse::parse_silent)) && - !(get_keyval(conf, "extended_v", vr, + !(get_keyval(conf, "extended_v", v_ext, 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"); } - x_reported = xr; + x_reported = x_ext; } else { x_reported = x; } @@ -1875,7 +1999,7 @@ std::istream & colvar::read_restart(std::istream &is) } if (is_enabled(f_cv_extended_Lagrangian)) { - v_reported = vr; + v_reported = v_ext; } else { v_reported = v_fdiff; } @@ -1901,8 +2025,8 @@ std::istream & colvar::read_traj(std::istream &is) } if (is_enabled(f_cv_extended_Lagrangian)) { - is >> xr; - x_reported = xr; + is >> x_ext; + x_reported = x_ext; } else { x_reported = x; } @@ -1913,8 +2037,8 @@ std::istream & colvar::read_traj(std::istream &is) is >> v_fdiff; if (is_enabled(f_cv_extended_Lagrangian)) { - is >> vr; - v_reported = vr; + is >> v_ext; + v_reported = v_ext; } else { v_reported = v_fdiff; } @@ -1955,11 +2079,11 @@ std::ostream & colvar::write_restart(std::ostream &os) { os << " extended_x " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << xr << "\n" + << x_ext << "\n" << " extended_v " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << vr << "\n"; + << v_ext << "\n"; } os << "}\n\n"; @@ -2190,6 +2314,7 @@ int colvar::calc_acf() acf_x_history_p = acf_x_history.begin(); break; + case acf_notset: default: break; } @@ -2222,6 +2347,7 @@ int colvar::calc_acf() history_incr(acf_x_history, acf_x_history_p); break; + case acf_notset: default: break; } @@ -2257,7 +2383,7 @@ void colvar::calc_vel_acf(std::list &v_list, void colvar::calc_coor_acf(std::list &x_list, - colvarvalue const &x) + colvarvalue const &x_now) { // same as above but for coordinates if (x_list.size() >= acf_length+acf_offset) { @@ -2269,7 +2395,7 @@ void colvar::calc_coor_acf(std::list &x_list, *(acf_i++) += x.norm2(); - colvarvalue::inner_opt(x, xs_i, x_list.end(), acf_i); + colvarvalue::inner_opt(x_now, xs_i, x_list.end(), acf_i); acf_nframes++; } @@ -2277,7 +2403,7 @@ void colvar::calc_coor_acf(std::list &x_list, void colvar::calc_p2coor_acf(std::list &x_list, - colvarvalue const &x) + colvarvalue const &x_now) { // same as above but with second order Legendre polynomial instead // of just the scalar product @@ -2291,7 +2417,7 @@ void colvar::calc_p2coor_acf(std::list &x_list, // value of P2(0) = 1 *(acf_i++) += 1.0; - colvarvalue::p2leg_opt(x, xs_i, x_list.end(), acf_i); + colvarvalue::p2leg_opt(x_now, xs_i, x_list.end(), acf_i); acf_nframes++; } @@ -2316,6 +2442,9 @@ int colvar::write_acf(std::ostream &os) case acf_p2coor: os << "Coordinate (2nd Legendre poly)"; break; + case acf_notset: + default: + break; } if (acf_colvar_name == name) { @@ -2420,7 +2549,7 @@ int colvar::calc_runave() << 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"; + << cvm::sqrt(runave_variance) << "\n"; } history_add_value(runave_length, *x_history_p, x); diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index a67749d577..74f7fdee51 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -92,7 +92,7 @@ public: static std::vector cv_features; /// \brief Implementation of the feature list accessor for colvar - virtual const std::vector &features() + virtual const std::vector &features() const { return cv_features; } @@ -133,7 +133,7 @@ protected: Here: S(x(t)) = x - s(t) = xr + s(t) = x_ext DS = Ds = delta */ @@ -170,13 +170,13 @@ protected: // Options for extended_lagrangian /// Restraint center - colvarvalue xr; + colvarvalue x_ext; /// Previous value of the restraint center; - colvarvalue prev_xr; + colvarvalue prev_x_ext; /// Velocity of the restraint center - colvarvalue vr; + colvarvalue v_ext; /// Previous velocity of the restraint center - colvarvalue prev_vr; + colvarvalue prev_v_ext; /// Mass of the restraint center cvm::real ext_mass; /// Restraint force constant @@ -273,6 +273,9 @@ public: /// Init output flags int init_output_flags(std::string const &conf); + /// \brief Initialize dependency tree + virtual int init_dependencies(); + private: /// Parse the CVC configuration for all components of a certain type template int init_components_type(std::string const &conf, @@ -373,7 +376,7 @@ protected: void update_active_cvc_square_norm(); /// \brief Absolute timestep number when this colvar was last updated - int prev_timestep; + cvm::step_number prev_timestep; public: @@ -383,32 +386,32 @@ public: /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags) inline size_t num_active_cvcs() const { return n_active_cvcs; } - /// \brief Use the internal metrics (as from \link cvc + /// \brief Use the internal metrics (as from \link colvar::cvc /// \endlink objects) to calculate square distances and gradients /// /// Handles correctly symmetries and periodic boundary conditions cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; - /// \brief Use the internal metrics (as from \link cvc + /// \brief Use the internal metrics (as from \link colvar::cvc /// \endlink objects) to calculate square distances and gradients /// /// Handles correctly symmetries and periodic boundary conditions colvarvalue dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const; - /// \brief Use the internal metrics (as from \link cvc + /// \brief Use the internal metrics (as from \link colvar::cvc /// \endlink objects) to calculate square distances and gradients /// /// Handles correctly symmetries and periodic boundary conditions colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; - /// \brief Use the internal metrics (as from \link cvc + /// \brief Use the internal metrics (as from \link colvar::cvc /// \endlink objects) to wrap a value into a standard interval /// /// Handles correctly symmetries and periodic boundary conditions - void wrap(colvarvalue &x) const; + void wrap(colvarvalue &x_unwrapped) const; /// Read the analysis tasks @@ -546,6 +549,7 @@ public: class polar_phi; class distance_inv; class distance_pairs; + class dipole_magnitude; class angle; class dipole_angle; class dihedral; @@ -574,7 +578,7 @@ public: protected: - /// \brief Array of \link cvc \endlink objects + /// \brief Array of \link colvar::cvc \endlink objects std::vector cvcs; /// \brief Flags to enable or disable cvcs at next colvar evaluation @@ -619,6 +623,9 @@ public: inline size_t n_components() const { return cvcs.size(); } + + /// \brief Get vector of vectors of atom IDs for all atom groups + virtual std::vector > get_atom_lists(); }; inline cvm::real const & colvar::force_constant() const @@ -655,6 +662,8 @@ inline colvarvalue const & colvar::total_force() const inline void colvar::add_bias_force(colvarvalue const &force) { + check_enabled(f_cv_gradient, + std::string("applying a force to the variable \""+name+"\"")); if (cvm::debug()) { cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n"); } diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h index 759b8d54a0..365f46148a 100644 --- a/lib/colvars/colvar_UIestimator.h +++ b/lib/colvars/colvar_UIestimator.h @@ -33,24 +33,24 @@ namespace UIestimator { public: n_matrix() {} - n_matrix(const std::vector & lowerboundary, // lowerboundary of x - const std::vector & upperboundary, // upperboundary of - const std::vector & width, // width of x - const int y_size) { // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered + n_matrix(const std::vector & lowerboundary_input, // lowerboundary of x + const std::vector & upperboundary_input, // upperboundary of + const std::vector & width_input, // width of x + const int y_size_input) { // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered int i; - this->lowerboundary = lowerboundary; - this->upperboundary = upperboundary; - this->width = width; - this->dimension = lowerboundary.size(); - this->y_size = y_size; // keep in mind the internal (spare) matrix is stored in diagonal form - this->y_total_size = int(std::pow(double(y_size), double(dimension)) + EPSILON); + this->lowerboundary = lowerboundary_input; + this->upperboundary = upperboundary_input; + this->width = width_input; + this->dimension = lowerboundary_input.size(); + this->y_size = y_size_input; // keep in mind the internal (spare) matrix is stored in diagonal form + this->y_total_size = int(cvm::pow(double(y_size_input), double(dimension)) + EPSILON); // the range of the matrix is [lowerboundary, upperboundary] x_total_size = 1; for (i = 0; i < dimension; i++) { - x_size.push_back(int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON)); + x_size.push_back(int((upperboundary_input[i] - lowerboundary_input[i]) / width_input[i] + EPSILON)); x_total_size *= x_size[i]; } @@ -89,9 +89,10 @@ namespace UIestimator { std::vector temp; // this vector is used in convert_x and convert_y to save computational resource - int i, j; - int convert_x(const std::vector & x) { // convert real x value to its interal index + + int i, j; + for (i = 0; i < dimension; i++) { temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON); } @@ -121,7 +122,7 @@ namespace UIestimator { int index = 0; for (i = 0; i < dimension; i++) { if (i + 1 < dimension) - index += temp[i] * int(std::pow(double(y_size), double(dimension - i - 1)) + EPSILON); + index += temp[i] * int(cvm::pow(double(y_size), double(dimension - i - 1)) + EPSILON); else index += temp[i]; } @@ -139,19 +140,19 @@ namespace UIestimator { public: n_vector() {} - n_vector(const std::vector & lowerboundary, // lowerboundary of x - const std::vector & upperboundary, // upperboundary of - const std::vector & width, // width of x - const int y_size, // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered + n_vector(const std::vector & lowerboundary_input, // lowerboundary of x + const std::vector & upperboundary_input, // upperboundary of + const std::vector & width_input, // width of x + const int y_size_input, // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered const T & default_value) { // the default value of T - this->width = width; - this->dimension = lowerboundary.size(); + this->width = width_input; + this->dimension = lowerboundary_input.size(); x_total_size = 1; for (int i = 0; i < dimension; i++) { - this->lowerboundary.push_back(lowerboundary[i] - (y_size - 1) / 2 * width[i] - EPSILON); - this->upperboundary.push_back(upperboundary[i] + (y_size - 1) / 2 * width[i] + EPSILON); + this->lowerboundary.push_back(lowerboundary_input[i] - (y_size_input - 1) / 2 * width_input[i] - EPSILON); + this->upperboundary.push_back(upperboundary_input[i] + (y_size_input - 1) / 2 * width_input[i] + EPSILON); x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + EPSILON)); x_total_size *= x_size[i]; @@ -215,26 +216,26 @@ namespace UIestimator { UIestimator() {} //called when (re)start an eabf simulation - UIestimator(const std::vector & lowerboundary, - const std::vector & upperboundary, - const std::vector & width, - const std::vector & krestr, // force constant in eABF - const std::string & output_filename, // the prefix of output files - const int output_freq, - const bool restart, // whether restart from a .count and a .grad file - const std::vector & input_filename, // the prefixes of input files - const double temperature) { + UIestimator(const std::vector & lowerboundary_input, + const std::vector & upperboundary_input, + const std::vector & width_input, + const std::vector & krestr_input, // force constant in eABF + const std::string & output_filename_input, // the prefix of output files + const int output_freq_input, + const bool restart_input, // whether restart from a .count and a .grad file + const std::vector & input_filename_input, // the prefixes of input files + const double temperature_input) { // initialize variables - this->lowerboundary = lowerboundary; - this->upperboundary = upperboundary; - this->width = width; - this->krestr = krestr; - this->output_filename = output_filename; - this->output_freq = output_freq; - this->restart = restart; - this->input_filename = input_filename; - this->temperature = temperature; + this->lowerboundary = lowerboundary_input; + this->upperboundary = upperboundary_input; + this->width = width_input; + this->krestr = krestr_input; + this->output_filename = output_filename_input; + this->output_freq = output_freq_input; + this->restart = restart_input; + this->input_filename = input_filename_input; + this->temperature = temperature_input; int i, j; @@ -300,7 +301,7 @@ namespace UIestimator { ~UIestimator() {} // called from MD engine every step - bool update(const int step, std::vector x, std::vector y) { + bool update(cvm::step_number step, std::vector x, std::vector y) { int i; @@ -431,7 +432,7 @@ namespace UIestimator { loop_flag_y[k] = loop_flag_x[k] - HALF_Y_SIZE * width[k]; } - int j = 0; + j = 0; while (j >= 0) { norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y); for (k = 0; k < dimension; k++) { @@ -672,7 +673,7 @@ namespace UIestimator { } // read input files - void read_inputfiles(const std::vector input_filename) + void read_inputfiles(const std::vector filename) { char sharp; double nothing; @@ -683,11 +684,11 @@ namespace UIestimator { std::vector position_temp(dimension, 0); std::vector grad_temp(dimension, 0); int count_temp = 0; - for (i = 0; i < int(input_filename.size()); i++) { + for (i = 0; i < int(filename.size()); i++) { int size = 1 , size_temp = 0; - std::string count_filename = input_filename[i] + ".UI.count"; - std::string grad_filename = input_filename[i] + ".UI.grad"; + std::string count_filename = filename[i] + ".UI.count"; + std::string grad_filename = filename[i] + ".UI.grad"; std::ifstream count_file(count_filename.c_str(), std::ios::in); std::ifstream grad_file(grad_filename.c_str(), std::ios::in); diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 3315007b54..eeb7985fec 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -22,7 +22,7 @@ cvm::atom::atom() index = -1; id = -1; mass = 1.0; - charge = 1.0; + charge = 0.0; reset_data(); } @@ -107,6 +107,8 @@ cvm::atom_group::~atom_group() delete fitting_group; fitting_group = NULL; } + + cvm::main()->unregister_named_atom_group(this); } @@ -183,10 +185,7 @@ int cvm::atom_group::init() // These may be overwritten by parse(), if a name is provided atoms.clear(); - - // TODO: check with proxy whether atom forces etc are available - init_ag_requires(); - + init_dependencies(); index = -1; b_dummy = false; @@ -207,8 +206,67 @@ int cvm::atom_group::init() } +int cvm::atom_group::init_dependencies() { + size_t i; + // Initialize static array once and for all + if (features().size() == 0) { + for (i = 0; i < f_ag_ntot; i++) { + modify_features().push_back(new feature); + } + + init_feature(f_ag_active, "active", f_type_dynamic); + init_feature(f_ag_center, "translational fit", f_type_static); + init_feature(f_ag_rotate, "rotational fit", f_type_static); + init_feature(f_ag_fitting_group, "fitting group", f_type_static); + init_feature(f_ag_explicit_gradient, "explicit atom gradient", f_type_dynamic); + init_feature(f_ag_fit_gradients, "fit gradients", f_type_user); + require_feature_self(f_ag_fit_gradients, f_ag_explicit_gradient); + + init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic); + + // parallel calculation implies that we have at least a scalable center of mass, + // but f_ag_scalable is kept as a separate feature to deal with future dependencies + init_feature(f_ag_scalable, "scalable group calculation", f_type_static); + init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static); + require_feature_self(f_ag_scalable, f_ag_scalable_com); + + // check that everything is initialized + for (i = 0; i < colvardeps::f_ag_ntot; i++) { + if (is_not_set(i)) { + cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description); + } + } + } + + // Initialize feature_states for each instance + // default as unavailable, not enabled + feature_states.reserve(f_ag_ntot); + for (i = 0; i < colvardeps::f_ag_ntot; i++) { + feature_states.push_back(feature_state(false, false)); + } + + // Features that are implemented (or not) by all atom groups + feature_states[f_ag_active].available = true; + // f_ag_scalable_com is provided by the CVC iff it is COM-based + feature_states[f_ag_scalable_com].available = false; + // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else) + feature_states[f_ag_scalable].available = true; + feature_states[f_ag_fit_gradients].available = true; + feature_states[f_ag_fitting_group].available = true; + feature_states[f_ag_explicit_gradient].available = true; + + return COLVARS_OK; +} + + int cvm::atom_group::setup() { + if (atoms_ids.size() == 0) { + atoms_ids.reserve(atoms.size()); + for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + atoms_ids.push_back(ai->id); + } + } for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { ai->update_mass(); ai->update_charge(); @@ -237,15 +295,6 @@ void cvm::atom_group::update_total_mass() } -void cvm::atom_group::reset_mass(std::string &name, int i, int j) -{ - update_total_mass(); - cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+ - cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+ - " atoms: total mass = "+cvm::to_str(total_mass)+".\n"); -} - - void cvm::atom_group::update_total_charge() { if (b_dummy) { @@ -264,6 +313,19 @@ void cvm::atom_group::update_total_charge() } +void cvm::atom_group::print_properties(std::string const &colvar_name, + int i, int j) +{ + if (cvm::proxy->updated_masses() && cvm::proxy->updated_charges()) { + cvm::log("Re-initialized atom group for variable \""+colvar_name+"\":"+ + cvm::to_str(i)+"/"+ + cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+ + " atoms: total mass = "+cvm::to_str(total_mass)+ + ", total charge = "+cvm::to_str(total_charge)+".\n"); + } +} + + int cvm::atom_group::parse(std::string const &group_conf) { cvm::log("Initializing atom group \""+key+"\".\n"); @@ -450,10 +512,21 @@ int cvm::atom_group::parse(std::string const &group_conf) if (cvm::debug()) cvm::log("Done initializing atom group \""+key+"\".\n"); - cvm::log("Atom group \""+key+"\" defined, "+ - cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+ - cvm::to_str(total_mass)+", total charge = "+ - cvm::to_str(total_charge)+".\n"); + { + std::string init_msg; + init_msg.append("Atom group \""+key+"\" defined with "+ + cvm::to_str(atoms_ids.size())+" atoms requested"); + if ((cvm::proxy)->updated_masses()) { + init_msg.append(": total mass = "+ + cvm::to_str(total_mass)); + if ((cvm::proxy)->updated_charges()) { + init_msg.append(", total charge = "+ + cvm::to_str(total_charge)); + } + } + init_msg.append(".\n"); + cvm::log(init_msg); + } if (b_print_atom_ids) { cvm::log("Internal definition of the atom group:\n"); @@ -464,7 +537,7 @@ int cvm::atom_group::parse(std::string const &group_conf) } -int cvm::atom_group::add_atoms_of_group(atom_group const * ag) +int cvm::atom_group::add_atoms_of_group(atom_group const *ag) { std::vector const &source_ids = ag->atoms_ids; @@ -696,6 +769,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf) return INPUT_ERROR; } } + enable(f_ag_fitting_group); } atom_group *group_for_fit = fitting_group ? fitting_group : this; @@ -800,24 +874,24 @@ int cvm::atom_group::create_sorted_ids() // Sort the internal IDs std::list sorted_atoms_ids_list; - for (size_t i = 0; i < this->size(); i++) { + for (size_t i = 0; i < atoms_ids.size(); i++) { sorted_atoms_ids_list.push_back(atoms_ids[i]); } sorted_atoms_ids_list.sort(); sorted_atoms_ids_list.unique(); - if (sorted_atoms_ids_list.size() != this->size()) { + if (sorted_atoms_ids_list.size() != atoms_ids.size()) { return cvm::error("Error: duplicate atom IDs in atom group? (found " + cvm::to_str(sorted_atoms_ids_list.size()) + " unique atom IDs instead of " + - cvm::to_str(this->size()) + ").\n", BUG_ERROR); + cvm::to_str(atoms_ids.size()) + ").\n", BUG_ERROR); } // Compute map between sorted and unsorted elements - sorted_atoms_ids.resize(this->size()); - sorted_atoms_ids_map.resize(this->size()); + sorted_atoms_ids.resize(atoms_ids.size()); + sorted_atoms_ids_map.resize(atoms_ids.size()); std::list::iterator lsii = sorted_atoms_ids_list.begin(); size_t ii = 0; - for ( ; ii < this->size(); lsii++, ii++) { + for ( ; ii < atoms_ids.size(); lsii++, ii++) { sorted_atoms_ids[ii] = *lsii; size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) - atoms_ids.begin(); @@ -1038,15 +1112,15 @@ int cvm::atom_group::calc_center_of_mass() } -int cvm::atom_group::calc_dipole(cvm::atom_pos const &com) +int cvm::atom_group::calc_dipole(cvm::atom_pos const &dipole_center) { if (b_dummy) { - cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR); - return COLVARS_ERROR; + return cvm::error("Error: trying to compute the dipole " + "of a dummy group.\n", INPUT_ERROR); } dip.reset(); for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { - dip += ai->charge * (ai->pos - com); + dip += ai->charge * (ai->pos - dipole_center); } return COLVARS_OK; } @@ -1056,13 +1130,12 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad) { if (b_dummy) return; - if (is_enabled(f_ag_scalable)) { - scalar_com_gradient = grad; - return; - } + scalar_com_gradient = grad; - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { - ai->grad = (ai->mass/total_mass) * grad; + if (!is_enabled(f_ag_scalable)) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + ai->grad = (ai->mass/total_mass) * grad; + } } } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 0b0dd62c70..9756e0e364 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -17,7 +17,7 @@ /// \brief Stores numeric id, mass and all mutable data for an atom, -/// mostly used by a \link cvc \endlink +/// mostly used by a \link colvar::cvc \endlink /// /// This class may be used to keep atomic data such as id, mass, /// position and collective variable derivatives) altogether. @@ -63,7 +63,7 @@ public: /// 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 + /// defined within the specific \link colvar::cvc \endlink /// implementation cvm::rvector grad; @@ -100,13 +100,19 @@ public: /// Get the latest value of the mass inline void update_mass() { - mass = (cvm::proxy)->get_atom_mass(index); + colvarproxy *p = cvm::proxy; + if (p->updated_masses()) { + mass = p->get_atom_mass(index); + } } /// Get the latest value of the charge inline void update_charge() { - charge = (cvm::proxy)->get_atom_charge(index); + colvarproxy *p = cvm::proxy; + if (p->updated_charges()) { + charge = p->get_atom_charge(index); + } } /// Get the current position @@ -145,7 +151,7 @@ public: /// \brief Group of \link atom \endlink objects, mostly used by a -/// \link cvc \endlink object to gather all atomic data +/// \link colvar::cvc \endlink object to gather all atomic data class colvarmodule::atom_group : public colvarparse, public colvardeps { @@ -174,6 +180,9 @@ public: /// \brief Set default values for common flags int init(); + /// \brief Initialize dependency tree + virtual int init_dependencies(); + /// \brief Update data required to calculate cvc's int setup(); @@ -198,16 +207,16 @@ public: /// \brief Remove an atom object from this group int remove_atom(cvm::atom_iter ai); - /// \brief Re-initialize the total mass of a group. + /// \brief Print the updated the total mass and charge of a group. /// This is needed in case the hosting MD code has an option to /// change atom masses after their initialization. - void reset_mass(std::string &name, int i, int j); + void print_properties(std::string const &colvar_name, int i, int j); /// \brief Implementation of the feature list for atom group static std::vector ag_features; /// \brief Implementation of the feature list accessor for atom group - virtual const std::vector &features() + virtual const std::vector &features() const { return ag_features; } @@ -347,15 +356,19 @@ public: /// Total mass of the atom group cvm::real total_mass; + + /// Update the total mass of the atom group void update_total_mass(); /// Total charge of the atom group cvm::real total_charge; + + /// Update the total mass of the group void update_total_charge(); /// \brief Don't apply any force on this group (use its coordinates /// only to calculate a colvar) - bool noforce; + bool noforce; /// \brief Get the current positions void read_positions(); @@ -423,20 +436,32 @@ public: /// \brief Calculate the center of mass of the atomic positions, assuming that /// they are already pbc-wrapped int calc_center_of_mass(); + private: + /// \brief Center of mass cvm::atom_pos com; + /// \brief The derivative of a scalar variable with respect to the COM // TODO for scalable calculations of more complex variables (e.g. rotation), // use a colvarvalue of vectors to hold the entire derivative cvm::rvector scalar_com_gradient; + public: - /// \brief Return the center of mass of the atomic positions + + /// \brief Return the center of mass (COM) of the atomic positions inline cvm::atom_pos center_of_mass() const { return com; } + /// \brief Return previously gradient of scalar variable with respect to the + /// COM + inline cvm::rvector center_of_mass_scalar_gradient() const + { + return scalar_com_gradient; + } + /// \brief Return a copy of the current atom positions, shifted by a constant vector std::vector positions_shifted(cvm::rvector const &shift) const; @@ -444,10 +469,15 @@ public: std::vector velocities() const; ///\brief Calculate the dipole of the atom group around the specified center - int calc_dipole(cvm::atom_pos const &com); + int calc_dipole(cvm::atom_pos const &dipole_center); + private: + + /// Dipole moment of the atom group cvm::rvector dip; + public: + ///\brief Return the (previously calculated) dipole of the atom group inline cvm::rvector dipole() const { diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 9363fcdcb6..724326d3b4 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -17,15 +17,14 @@ colvarbias::colvarbias(char const *key) : bias_type(to_lower_cppstr(key)) { - init_cvb_requires(); - + description = "uninitialized " + cvm::to_str(key) + " bias"; + init_dependencies(); rank = 1; has_data = false; b_output_energy = false; reset(); - state_file_step = 0; - description = "uninitialized " + cvm::to_str(key) + " bias"; + state_file_step = 0L; } @@ -76,6 +75,7 @@ int colvarbias::init(std::string const &conf) cvm::error("Error: no collective variables specified.\n", INPUT_ERROR); return INPUT_ERROR; } + } else { cvm::log("Reinitializing bias \""+name+"\".\n"); } @@ -98,6 +98,70 @@ int colvarbias::init(std::string const &conf) } +int colvarbias::init_dependencies() { + int i; + if (features().size() == 0) { + for (i = 0; i < f_cvb_ntot; i++) { + modify_features().push_back(new feature); + } + + init_feature(f_cvb_active, "active", f_type_dynamic); + require_feature_children(f_cvb_active, f_cv_active); + + init_feature(f_cvb_awake, "awake", f_type_static); + require_feature_self(f_cvb_awake, f_cvb_active); + + init_feature(f_cvb_apply_force, "apply force", f_type_user); + require_feature_children(f_cvb_apply_force, f_cv_gradient); + + init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic); + require_feature_children(f_cvb_get_total_force, f_cv_total_force); + + init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user); + require_feature_self(f_cvb_output_acc_work, f_cvb_apply_force); + + init_feature(f_cvb_history_dependent, "history-dependent", f_type_static); + + init_feature(f_cvb_time_dependent, "time-dependent", f_type_static); + + init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static); + require_feature_children(f_cvb_scalar_variables, f_cv_scalar); + + init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static); + + init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic); + require_feature_self(f_cvb_calc_ti_samples, f_cvb_get_total_force); + require_feature_children(f_cvb_calc_ti_samples, f_cv_grid); + + init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user); + require_feature_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples); + + init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user); + require_feature_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples); + + // check that everything is initialized + for (i = 0; i < colvardeps::f_cvb_ntot; i++) { + if (is_not_set(i)) { + cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description); + } + } + } + + // Initialize feature_states for each instance + feature_states.reserve(f_cvb_ntot); + for (i = 0; i < f_cvb_ntot; i++) { + feature_states.push_back(feature_state(true, false)); + // Most features are available, so we set them so + // and list exceptions below + } + + // only compute TI samples when deriving from colvarbias_ti + feature_states[f_cvb_calc_ti_samples].available = false; + + return COLVARS_OK; +} + + int colvarbias::reset() { bias_energy = 0.0; @@ -217,6 +281,9 @@ int colvarbias::update() void colvarbias::communicate_forces() { + if (! is_enabled(f_cvb_apply_force)) { + return; + } size_t i = 0; for (i = 0; i < num_variables(); i++) { if (cvm::debug()) { @@ -345,7 +412,8 @@ std::istream & colvarbias::read_state(std::istream &is) (set_state_params(conf) != COLVARS_OK) ) { cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ - cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n", INPUT_ERROR); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -355,7 +423,8 @@ std::istream & colvarbias::read_state(std::istream &is) if (!read_state_data(is)) { cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ - cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n", INPUT_ERROR); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -365,7 +434,8 @@ std::istream & colvarbias::read_state(std::istream &is) if (brace != "}") { cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+ this->name+"\": no matching brace at position "+ - cvm::to_str(is.tellg())+" in stream.\n"); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n"); is.setstate(std::ios::failbit); } @@ -381,7 +451,8 @@ std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key !(key_in == to_lower_cppstr(std::string(key))) ) { cvm::error("Error: in reading restart configuration for "+ bias_type+" bias \""+this->name+"\" at position "+ - cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n", INPUT_ERROR); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -640,7 +711,7 @@ int colvarbias_ti::write_output_files() cvm::proxy->close_output_stream(ti_count_file_name); } - std::string const ti_grad_file_name(ti_output_prefix+".ti.grad"); + std::string const ti_grad_file_name(ti_output_prefix+".ti.force"); os = cvm::proxy->output_stream(ti_grad_file_name); if (os) { ti_avg_forces->write_multicol(*os); diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 391826e79e..eac88a7f18 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -96,6 +96,9 @@ public: /// \brief Parse config string and (re)initialize virtual int init(std::string const &conf); + /// \brief Initialize dependency tree + virtual int init_dependencies(); + /// \brief Set to zero all mutable data virtual int reset(); @@ -181,7 +184,7 @@ public: static std::vector cvb_features; /// \brief Implementation of the feature list accessor for colvarbias - virtual const std::vector &features() + virtual const std::vector &features() const { return cvb_features; } @@ -220,7 +223,7 @@ protected: bool has_data; /// \brief Step number read from the last state file - size_t state_file_step; + cvm::step_number state_file_step; }; diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index e5edd92ae6..bac021be99 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -191,10 +191,10 @@ int colvarbias_abf::init(std::string const &conf) // Projected ABF get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0); // Parameters for integrating initial (and final) gradient data - get_keyval(conf, "integrateInitSteps", integrate_initial_steps, 1e4); + get_keyval(conf, "integrateInitMaxIterations", integrate_initial_iterations, 1e4); get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6); // for updating the integrated PMF on the fly - get_keyval(conf, "integrateSteps", integrate_steps, 100); + get_keyval(conf, "integrateMaxIterations", integrate_iterations, 100); get_keyval(conf, "integrateTol", integrate_tol, 1e-4); } } else { @@ -366,10 +366,10 @@ int colvarbias_abf::update() if ( b_integrate ) { if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) { cvm::real err; - int iter = pmf->integrate(integrate_steps, integrate_tol, err); - if ( iter == integrate_steps ) { + int iter = pmf->integrate(integrate_iterations, integrate_tol, err); + if ( iter == integrate_iterations ) { cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol) - + " in " + cvm::to_str(integrate_steps) + + " in " + cvm::to_str(integrate_iterations) + " steps. Residual error: " + cvm::to_str(err)); } pmf->set_zero_minimum(); // TODO: do this only when necessary @@ -597,7 +597,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app if (b_integrate) { // Do numerical integration (to high precision) and output a PMF cvm::real err; - pmf->integrate(integrate_initial_steps, integrate_initial_tol, err); + pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err); pmf->set_zero_minimum(); std::string pmf_out_name = prefix + ".pmf"; @@ -661,7 +661,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app // Do numerical integration (to high precision) and output a PMF cvm::real err; czar_pmf->set_div(); - czar_pmf->integrate(integrate_initial_steps, integrate_initial_tol, err); + czar_pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err); czar_pmf->set_zero_minimum(); std::string czar_pmf_out_name = prefix + ".czar.pmf"; diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index 52bf2df210..4bcc149da5 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -27,9 +27,13 @@ class colvarbias_abf : public colvarbias { public: + /// Constructor for ABF bias colvarbias_abf(char const *key); + /// Initializer for ABF bias virtual int init(std::string const &conf); + /// Default destructor for ABF bias virtual ~colvarbias_abf(); + /// Per-timestep update of ABF bias virtual int update(); private: @@ -40,11 +44,17 @@ private: /// Base filename(s) for reading previous gradient data (replaces data from restart file) std::vector input_prefix; + /// Adapt the bias at each time step (as opposed to keeping it constant)? bool update_bias; + /// Use normalized definition of PMF for distance functions? (flat at long distances) + /// by including the Jacobian term separately of the recorded PMF bool hide_Jacobian; + /// Integrate gradients into a PMF on output bool b_integrate; + /// Number of samples per bin before applying the full biasing force size_t full_samples; + /// Number of samples per bin before applying a scaled-down biasing force size_t min_samples; /// frequency for updating output files int output_freq; @@ -52,6 +62,7 @@ private: bool b_history_files; /// Write CZAR output file for stratified eABF (.zgrad) bool b_czar_window_file; + /// Number of timesteps between recording data in history files (if non-zero) size_t history_freq; /// Umbrella Integration estimator of free energy from eABF UIestimator::UIestimator eabf_UI; @@ -63,25 +74,30 @@ private: /// Frequency for updating pABF PMF (if zero, pABF is not used) int pabf_freq; /// Max number of CG iterations for integrating PMF at startup and for file output - int integrate_initial_steps; + int integrate_initial_iterations; /// Tolerance for integrating PMF at startup and for file output cvm::real integrate_initial_tol; /// Max number of CG iterations for integrating PMF at on-the-fly pABF updates - int integrate_steps; + int integrate_iterations; /// Tolerance for integrating PMF at on-the-fly pABF updates cvm::real integrate_tol; - /// Cap the biasing force to be applied? + /// Cap the biasing force to be applied? (option maxForce) bool cap_force; + /// Maximum force to be applied std::vector max_force; - // Frequency for updating 2D gradients - int integrate_freq; - // Internal data and methods - std::vector bin, force_bin, z_bin; - gradient_t system_force, applied_force; + /// Current bin in sample grid + std::vector bin; + /// Current bin in force grid + std::vector force_bin; + /// Cuurent bin in "actual" coordinate, when running extended Lagrangian dynamics + std::vector z_bin; + + /// Measured instantaneous system force + gradient_t system_force; /// n-dim grid of free energy gradients colvar_grid_gradient *gradients; @@ -118,7 +134,7 @@ private: // shared ABF bool shared_on; size_t shared_freq; - int shared_last_step; + cvm::step_number shared_last_step; // Share between replicas -- may be called independently of update virtual int replica_share(); diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h index b9c1b49950..23565caa5c 100644 --- a/lib/colvars/colvarbias_histogram.h +++ b/lib/colvars/colvarbias_histogram.h @@ -37,7 +37,7 @@ protected: std::string out_name, out_name_dx; size_t output_freq; - /// If one or more of the variables are \link type_vector \endlink, treat them as arrays of this length + /// If one or more of the variables are \link colvarvalue::type_vector \endlink, treat them as arrays of this length size_t colvar_array_size; /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram std::vector weights; diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index f3ae3631a0..27781ec733 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -11,7 +11,6 @@ #include #include #include -#include #include // used to set the absolute path of a replica file @@ -39,6 +38,8 @@ colvarbias_meta::colvarbias_meta(char const *key) new_hills_begin = hills.end(); hills_traj_os = NULL; replica_hills_os = NULL; + + ebmeta_equil_steps = 0L; } @@ -61,7 +62,7 @@ int colvarbias_meta::init(std::string const &conf) enable(f_cvb_history_dependent); } - get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0); + get_keyval(conf, "hillWidth", hill_width, cvm::sqrt(2.0 * PI) / 2.0); cvm::log("Half-widths of the Gaussian hills (sigma's):\n"); for (size_t i = 0; i < num_variables(); i++) { cvm::log(variables(i)->name+std::string(": ")+ @@ -201,6 +202,7 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf) } target_dist = new colvar_grid_scalar(); target_dist->init_from_colvars(colvars); + std::string target_dist_file; get_keyval(conf, "targetdistfile", target_dist_file); std::ifstream targetdiststream(target_dist_file.c_str()); target_dist->read_multicol(targetdiststream); @@ -221,9 +223,9 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf) } // normalize target distribution and multiply by effective volume = exp(differential entropy) target_dist->multiply_constant(1.0/target_dist->integral()); - cvm::real volume = std::exp(target_dist->entropy()); + cvm::real volume = cvm::exp(target_dist->entropy()); target_dist->multiply_constant(volume); - get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0); + get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, ebmeta_equil_steps); } return COLVARS_OK; @@ -291,7 +293,7 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h) // 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) { + if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) { hills_off_grid.push_back(h); } } @@ -387,7 +389,7 @@ int colvarbias_meta::update_grid_params() // first of all, expand the grids, if specified bool changed_grids = false; int const min_buffer = - (3 * (size_t) std::floor(hill_width)) + 1; + (3 * (size_t) cvm::floor(hill_width)) + 1; std::vector new_sizes(hills_energy->sizes()); std::vector new_lower_boundaries(hills_energy->lower_boundaries); @@ -492,9 +494,9 @@ int colvarbias_meta::update_bias() if (ebmeta) { hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index()); - if(cvm::step_absolute() <= long(ebmeta_equil_steps)) { + if(cvm::step_absolute() <= ebmeta_equil_steps) { cvm::real const hills_lambda = - (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) / + (cvm::real(ebmeta_equil_steps - cvm::step_absolute())) / (cvm::real(ebmeta_equil_steps)); hills_scale = hills_lambda + (1-hills_lambda)*hills_scale; } @@ -508,7 +510,7 @@ int colvarbias_meta::update_bias() } else { calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here); } - hills_scale *= std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); + hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); } switch (comm) { @@ -710,7 +712,7 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, // 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)); + h->value(cvm::exp(-0.5*cv_sqdev)); } energy += h->energy(); } @@ -904,7 +906,7 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first 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) { + if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) { hills_off_grid.push_back(*h); } } @@ -1427,8 +1429,8 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) return is; } - size_t h_it; - get_keyval(data, "step", h_it, 0, parse_silent); + cvm::step_number h_it; + get_keyval(data, "step", h_it, 0L, parse_silent); if (h_it <= state_file_step) { if (cvm::debug()) cvm::log("Skipping a hill older than the state file for metadynamics bias \""+ @@ -1457,7 +1459,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) std::vector h_widths(num_variables()); get_keyval(data, "widths", h_widths, - std::vector(num_variables(), (std::sqrt(2.0 * PI) / 2.0)), + std::vector(num_variables(), (cvm::sqrt(2.0 * PI) / 2.0)), parse_silent); std::string h_replica = ""; @@ -1482,7 +1484,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) // 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) { + if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) { hills_off_grid.push_back(hills.back()); } } diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 78b2d35d41..0ba2bef1c3 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -19,8 +19,8 @@ #include "colvargrid.h" /// Metadynamics bias (implementation of \link colvarbias \endlink) -class colvarbias_meta - : public virtual colvarbias, +class colvarbias_meta + : public virtual colvarbias, public virtual colvarbias_ti { @@ -174,12 +174,14 @@ protected: /// \brief Biasing temperature in well-tempered metadynamics cvm::real bias_temperature; - // EBmeta parameters + /// Ensemble-biased metadynamics (EBmeta) flag bool ebmeta; + + /// Target distribution for EBmeta colvar_grid_scalar* target_dist; - std::string target_dist_file; - cvm::real target_dist_volume; - size_t ebmeta_equil_steps; + + /// Number of equilibration steps for EBmeta + cvm::step_number ebmeta_equil_steps; /// \brief Try to read the restart information by allocating new @@ -285,7 +287,7 @@ public: friend class colvarbias_meta; /// Time step at which this hill was added - size_t it; + cvm::step_number it; /// Identity of the replica who added this hill (only in multiple replica simulations) std::string replica; @@ -296,9 +298,9 @@ public: /// replica (optional) Identity of the replica which creates the /// hill inline hill(cvm::real const &W_in, - std::vector &cv, - cvm::real const &hill_width, - std::string const &replica_in = "") + std::vector &cv, + cvm::real const &hill_width, + std::string const &replica_in = "") : sW(1.0), W(W_in), centers(cv.size()), @@ -325,11 +327,11 @@ public: /// weight Weight of the hill \param centers Center of the hill /// \param widths Width of the hill around centers \param replica /// (optional) Identity of the replica which creates the hill - inline hill(size_t const &it_in, - cvm::real const &W_in, - std::vector const ¢ers_in, - std::vector const &widths_in, - std::string const &replica_in = "") + inline hill(cvm::step_number const &it_in, + cvm::real const &W_in, + std::vector const ¢ers_in, + std::vector const &widths_in, + std::string const &replica_in = "") : sW(1.0), W(W_in), centers(centers_in), diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 2daf7a0876..90588f5a1f 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -7,8 +7,6 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include - #include "colvarmodule.h" #include "colvarproxy.h" #include "colvarvalue.h" @@ -179,7 +177,7 @@ int colvarbias_restraint_k::change_configuration(std::string const &conf) colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key) { target_nstages = 0; - target_nsteps = 0; + target_nsteps = 0L; stage = 0; acc_work = 0.0; b_chg_centers = false; @@ -241,8 +239,8 @@ int colvarbias_restraint_moving::set_state_params(std::string const &conf) { if (b_chg_centers || b_chg_force_k) { if (target_nstages) { - if (!get_keyval(conf, "stage", stage)) - cvm::error("Error: current stage is missing from the restart.\n"); + get_keyval(conf, "stage", stage, stage, + colvarparse::parse_restart | colvarparse::parse_required); } } return COLVARS_OK; @@ -436,13 +434,13 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con colvarbias_restraint::set_state_params(conf); if (b_chg_centers) { - // cvm::log ("Reading the updated restraint centers from the restart.\n"); - if (!get_keyval(conf, "centers", colvar_centers)) - cvm::error("Error: restraint centers are missing from the restart.\n"); - if (is_enabled(f_cvb_output_acc_work)) { - if (!get_keyval(conf, "accumulatedWork", acc_work)) - cvm::error("Error: accumulatedWork is missing from the restart.\n"); - } + get_keyval(conf, "centers", colvar_centers, colvar_centers, + colvarparse::parse_restart | colvarparse::parse_required); + } + + if (is_enabled(f_cvb_output_acc_work)) { + get_keyval(conf, "accumulatedWork", acc_work, acc_work, + colvarparse::parse_restart | colvarparse::parse_required); } return COLVARS_OK; @@ -563,7 +561,7 @@ int colvarbias_restraint_k_moving::update() lambda = 0.0; } force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow(lambda, force_k_exp); + * cvm::pow(lambda, force_k_exp); cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda) + ", k = " + cvm::to_str(force_k)); @@ -585,7 +583,7 @@ int colvarbias_restraint_k_moving::update() for (size_t i = 0; i < num_variables(); i++) { dU_dk += d_restraint_potential_dk(i); } - restraint_FE += force_k_exp * std::pow(lambda, force_k_exp - 1.0) + restraint_FE += force_k_exp * cvm::pow(lambda, force_k_exp - 1.0) * (target_force_k - starting_force_k) * dU_dk; } @@ -608,7 +606,7 @@ int colvarbias_restraint_k_moving::update() lambda = cvm::real(stage) / cvm::real(target_nstages); } force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow(lambda, force_k_exp); + * cvm::pow(lambda, force_k_exp); cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda) + ", k = " + cvm::to_str(force_k)); @@ -622,7 +620,7 @@ int colvarbias_restraint_k_moving::update() lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps); cvm::real const force_k_old = force_k; force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow(lambda, force_k_exp); + * cvm::pow(lambda, force_k_exp); force_k_incr = force_k - force_k_old; } } @@ -672,13 +670,13 @@ int colvarbias_restraint_k_moving::set_state_params(std::string const &conf) colvarbias_restraint::set_state_params(conf); if (b_chg_force_k) { - // cvm::log ("Reading the updated force constant from the restart.\n"); - if (!get_keyval(conf, "forceConstant", force_k, force_k)) - cvm::error("Error: force constant is missing from the restart.\n"); - if (is_enabled(f_cvb_output_acc_work)) { - if (!get_keyval(conf, "accumulatedWork", acc_work)) - cvm::error("Error: accumulatedWork is missing from the restart.\n"); - } + get_keyval(conf, "forceConstant", force_k, force_k, + colvarparse::parse_restart | colvarparse::parse_required); + } + + if (is_enabled(f_cvb_output_acc_work)) { + get_keyval(conf, "accumulatedWork", acc_work, acc_work, + colvarparse::parse_restart | colvarparse::parse_required); } return COLVARS_OK; @@ -719,7 +717,8 @@ std::istream & colvarbias_restraint::read_state(std::istream &is) (set_state_params(conf) != COLVARS_OK) ) { cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ - cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n", INPUT_ERROR); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -729,7 +728,8 @@ std::istream & colvarbias_restraint::read_state(std::istream &is) if (!read_state_data(is)) { cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ - cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); + cvm::to_str(static_cast(is.tellg()))+ + " in stream.\n", INPUT_ERROR); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -740,7 +740,7 @@ std::istream & colvarbias_restraint::read_state(std::istream &is) cvm::log("brace = "+brace+"\n"); cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+ this->name+"\": no matching brace at position "+ - cvm::to_str(is.tellg())+" in stream.\n"); + cvm::to_str(static_cast(is.tellg()))+" in stream.\n"); is.setstate(std::ios::failbit); } @@ -787,11 +787,11 @@ int colvarbias_restraint_harmonic::init(std::string const &conf) colvarbias_restraint_k_moving::init(conf); for (size_t i = 0; i < num_variables(); i++) { - if (variables(i)->width != 1.0) - cvm::log("The force constant for colvar \""+variables(i)->name+ - "\" will be rescaled to "+ - cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+ - " according to the specified width.\n"); + cvm::real const w = variables(i)->width; + cvm::log("The force constant for colvar \""+variables(i)->name+ + "\" will be rescaled to "+ + cvm::to_str(force_k/(w*w))+ + " according to the specified width ("+cvm::to_str(w)+").\n"); } return COLVARS_OK; @@ -1014,7 +1014,7 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) INPUT_ERROR); return INPUT_ERROR; } - force_k = std::sqrt(lower_wall_k * upper_wall_k); + force_k = cvm::sqrt(lower_wall_k * upper_wall_k); // transform the two constants to relative values using gemetric mean as ref // to preserve force_k if provided as single parameter // (allow changing both via force_k) @@ -1037,25 +1037,21 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) if (lower_walls.size() > 0) { for (i = 0; i < num_variables(); i++) { - if (variables(i)->width != 1.0) - cvm::log("The lower wall force constant for colvar \""+ - variables(i)->name+ - "\" will be rescaled to "+ - cvm::to_str(lower_wall_k * force_k / - (variables(i)->width * variables(i)->width))+ - " according to the specified width.\n"); + cvm::real const w = variables(i)->width; + cvm::log("The lower wall force constant for colvar \""+ + variables(i)->name+"\" will be rescaled to "+ + cvm::to_str(lower_wall_k * force_k / (w*w))+ + " according to the specified width ("+cvm::to_str(w)+").\n"); } } if (upper_walls.size() > 0) { for (i = 0; i < num_variables(); i++) { - if (variables(i)->width != 1.0) - cvm::log("The upper wall force constant for colvar \""+ - variables(i)->name+ - "\" will be rescaled to "+ - cvm::to_str(upper_wall_k * force_k / - (variables(i)->width * variables(i)->width))+ - " according to the specified width.\n"); + cvm::real const w = variables(i)->width; + cvm::log("The upper wall force constant for colvar \""+ + variables(i)->name+"\" will be rescaled to "+ + cvm::to_str(upper_wall_k * force_k / (w*w))+ + " according to the specified width ("+cvm::to_str(w)+").\n"); } } @@ -1225,11 +1221,11 @@ int colvarbias_restraint_linear::init(std::string const &conf) INPUT_ERROR); return INPUT_ERROR; } - if (variables(i)->width != 1.0) - cvm::log("The force constant for colvar \""+variables(i)->name+ - "\" will be rescaled to "+ - cvm::to_str(force_k / variables(i)->width)+ - " according to the specified width.\n"); + cvm::real const w = variables(i)->width; + cvm::log("The force constant for colvar \""+variables(i)->name+ + "\" will be rescaled to "+ + cvm::to_str(force_k / w)+ + " according to the specified width ("+cvm::to_str(w)+").\n"); } return COLVARS_OK; @@ -1367,6 +1363,7 @@ colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key) int colvarbias_restraint_histogram::init(std::string const &conf) { colvarbias::init(conf); + enable(f_cvb_apply_force); get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary); get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary); @@ -1390,7 +1387,7 @@ int colvarbias_restraint_histogram::init(std::string const &conf) cvm::real const nbins = (upper_boundary - lower_boundary) / width; int const nbins_round = (int)(nbins); - if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) { + if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) { cvm::log("Warning: grid interval ("+ cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+ cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+ @@ -1440,7 +1437,7 @@ int colvarbias_restraint_histogram::init(std::string const &conf) } } cvm::real const ref_integral = ref_p.sum() * width; - if (std::fabs(ref_integral - 1.0) > 1.0e-03) { + if (cvm::fabs(ref_integral - 1.0) > 1.0e-03) { cvm::log("Reference distribution not normalized, normalizing to unity.\n"); ref_p /= ref_integral; } @@ -1471,7 +1468,7 @@ int colvarbias_restraint_histogram::update() vector_size += variables(icv)->value().size(); } - cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size); + cvm::real const norm = 1.0/(cvm::sqrt(2.0*PI)*gaussian_width*vector_size); // calculate the histogram p.reset(); @@ -1482,7 +1479,7 @@ int colvarbias_restraint_histogram::update() size_t igrid; for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width); - p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / + p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / (2.0 * gaussian_width * gaussian_width)); } } else if (cv.type() == colvarvalue::type_vector) { @@ -1492,7 +1489,7 @@ int colvarbias_restraint_histogram::update() size_t igrid; for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width); - p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / + p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / (2.0 * gaussian_width * gaussian_width)); } } @@ -1523,7 +1520,7 @@ int colvarbias_restraint_histogram::update() for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width); force += force_k_cv * p_diff[igrid] * - norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / + norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / (2.0 * gaussian_width * gaussian_width)) * (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width)); } @@ -1536,7 +1533,7 @@ int colvarbias_restraint_histogram::update() for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width); force += force_k_cv * p_diff[igrid] * - norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / + norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) / (2.0 * gaussian_width * gaussian_width)) * (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width)); } @@ -1550,7 +1547,7 @@ int colvarbias_restraint_histogram::update() } -std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) +int colvarbias_restraint_histogram::write_output_files() { if (b_write_histogram) { std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat"); @@ -1558,6 +1555,9 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) *os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width) << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3) << ")\n"; + + os->setf(std::ios::fixed, std::ios::floatfield); + size_t igrid; for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+1)*width); @@ -1572,13 +1572,7 @@ std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) } cvm::proxy->close_output_stream(file_name); } - return os; -} - - -std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is) -{ - return is; + return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index 3ee999c262..6493f7f16b 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -132,7 +132,7 @@ protected: /// \brief Number of steps required to reach the target force constant /// or restraint centers - long target_nsteps; + cvm::step_number target_nsteps; /// \brief Accumulated work (computed when outputAccumulatedWork == true) cvm::real acc_work; @@ -328,8 +328,7 @@ public: virtual int update(); - virtual std::istream & read_restart(std::istream &is); - virtual std::ostream & write_restart(std::ostream &os); + virtual int write_output_files(); virtual std::ostream & write_traj_label(std::ostream &os); virtual std::ostream & write_traj(std::ostream &os); diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index cb272eed05..3075ed82ca 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -20,7 +20,8 @@ colvar::cvc::cvc() b_periodic(false), b_try_scalable(true) { - init_cvc_requires(); + description = "uninitialized colvar component"; + init_dependencies(); sup_coeff = 1.0; period = 0.0; wrap_center = 0.0; @@ -33,7 +34,8 @@ colvar::cvc::cvc(std::string const &conf) b_periodic(false), b_try_scalable(true) { - init_cvc_requires(); + description = "uninitialized colvar component"; + init_dependencies(); sup_coeff = 1.0; period = 0.0; wrap_center = 0.0; @@ -176,6 +178,100 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, } +int colvar::cvc::init_dependencies() { + size_t i; + // Initialize static array once and for all + if (features().size() == 0) { + for (i = 0; i < colvardeps::f_cvc_ntot; i++) { + modify_features().push_back(new feature); + } + + init_feature(f_cvc_active, "active", f_type_dynamic); +// The dependency below may become useful if we use dynamic atom groups +// require_feature_children(f_cvc_active, f_ag_active); + + init_feature(f_cvc_scalar, "scalar", f_type_static); + + init_feature(f_cvc_gradient, "gradient", f_type_dynamic); + + init_feature(f_cvc_explicit_gradient, "explicit gradient", f_type_static); + require_feature_children(f_cvc_explicit_gradient, f_ag_explicit_gradient); + + init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic); + require_feature_self(f_cvc_inv_gradient, f_cvc_gradient); + + init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user); + require_feature_self(f_cvc_debug_gradient, f_cvc_gradient); + require_feature_self(f_cvc_debug_gradient, f_cvc_explicit_gradient); + + init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic); + require_feature_self(f_cvc_Jacobian, f_cvc_inv_gradient); + + init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static); + + init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user); + + // Compute total force on first site only to avoid unwanted + // coupling to other colvars (see e.g. Ciccotti et al., 2005) + init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user); + require_feature_self(f_cvc_one_site_total_force, f_cvc_com_based); + + init_feature(f_cvc_scalable, "scalable calculation", f_type_static); + require_feature_self(f_cvc_scalable, f_cvc_scalable_com); + + init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static); + require_feature_self(f_cvc_scalable_com, f_cvc_com_based); + + + // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group + // require_feature_children(f_cvc_scalable, f_ag_scalable); + // require_feature_children(f_cvc_scalable_com, f_ag_scalable_com); + + // check that everything is initialized + for (i = 0; i < colvardeps::f_cvc_ntot; i++) { + if (is_not_set(i)) { + cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description); + } + } + } + + // Initialize feature_states for each instance + // default as available, not enabled + // except dynamic features which default as unavailable + feature_states.reserve(f_cvc_ntot); + for (i = 0; i < colvardeps::f_cvc_ntot; i++) { + bool avail = is_dynamic(i) ? false : true; + feature_states.push_back(feature_state(avail, false)); + } + + // Features that are implemented by all cvcs by default + // Each cvc specifies what other features are available + feature_states[f_cvc_active].available = true; + feature_states[f_cvc_gradient].available = true; + + // CVCs are enabled from the start - get disabled based on flags + enable(f_cvc_active); + // feature_states[f_cvc_active].enabled = true; + + // Explicit gradients are implemented in mosts CVCs. Exceptions must be specified explicitly. + // feature_states[f_cvc_explicit_gradient].enabled = true; + enable(f_cvc_explicit_gradient); + + // Use minimum-image distances by default + // feature_states[f_cvc_pbc_minimum_image].enabled = true; + enable(f_cvc_pbc_minimum_image); + + // Features that are implemented by default if their requirements are + feature_states[f_cvc_one_site_total_force].available = true; + + // Features That are implemented only for certain simulation engine configurations + feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK); + feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available; + + return COLVARS_OK; +} + + int colvar::cvc::setup() { description = "cvc " + name; @@ -192,6 +288,7 @@ colvar::cvc::~cvc() } } + void colvar::cvc::read_data() { size_t ig; @@ -214,6 +311,66 @@ void colvar::cvc::read_data() } +std::vector > colvar::cvc::get_atom_lists() +{ + std::vector > lists; + + std::vector::iterator agi = atom_groups.begin(); + for ( ; agi != atom_groups.end(); ++agi) { + (*agi)->create_sorted_ids(); + lists.push_back((*agi)->sorted_ids()); + if ((*agi)->is_enabled(f_ag_fitting_group) && (*agi)->is_enabled(f_ag_fit_gradients)) { + cvm::atom_group &fg = *((*agi)->fitting_group); + fg.create_sorted_ids(); + lists.push_back(fg.sorted_ids()); + } + } + return lists; +} + + +void colvar::cvc::collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients) +{ + // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx + cvm::real coeff = sup_coeff * cvm::real(sup_np) * + cvm::integer_power(value().real_value, sup_np-1); + + for (size_t j = 0; j < atom_groups.size(); j++) { + + cvm::atom_group &ag = *(atom_groups[j]); + + // If necessary, apply inverse rotation to get atomic + // gradient in the laboratory frame + if (ag.b_rotate) { + cvm::rotation const rot_inv = ag.rot.inverse(); + + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad); + } + + } else { + + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; + } + } + if (ag.is_enabled(f_ag_fitting_group) && ag.is_enabled(f_ag_fit_gradients)) { + cvm::atom_group const &fg = *(ag.fitting_group); + for (size_t k = 0; k < fg.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + fg[k].id) - atom_ids.begin(); + // fit gradients are in the unrotated (simulation) frame + atomic_gradients[a] += coeff * fg.fit_gradients[k]; + } + } + } +} + + void colvar::cvc::calc_force_invgrads() { cvm::error("Error: calculation of inverse gradients is not implemented " @@ -306,8 +463,8 @@ void colvar::cvc::debug_gradients() 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"); + cvm::to_str(cvm::fabs(x_1 - x_0 - dx_pred) / + cvm::fabs(x_1 - x_0), 12, 5)+"\n"); } } @@ -341,8 +498,8 @@ void colvar::cvc::debug_gradients() 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), + cvm::to_str(cvm::fabs (x_1 - x_0 - dx_pred) / + cvm::fabs (x_1 - x_0), 12, 5)+ ".\n"); } @@ -378,7 +535,7 @@ colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1, } -void colvar::cvc::wrap(colvarvalue &x) const +void colvar::cvc::wrap(colvarvalue &x_unwrapped) const { return; } diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 1a6df3771e..f615680ba6 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -27,12 +27,12 @@ /// \brief Colvar component (base class for collective variables) /// -/// A \link cvc \endlink object (or an object of a +/// A \link colvar::cvc \endlink object (or an object of a /// cvc-derived class) implements the calculation of a collective /// variable, its gradients and any other related physical quantities /// that depend on microscopic degrees of freedom. /// -/// No restriction is set to what kind of calculation a \link cvc \endlink +/// No restriction is set to what kind of calculation a \link colvar::cvc \endlink /// object performs (usually an analytical function of atomic coordinates). /// The only constraints are that: \par /// @@ -42,9 +42,9 @@ /// alike, and allows an automatic selection of the applicable algorithms. /// /// - The object provides an implementation \link apply_force() \endlink to -/// apply forces to atoms. Typically, one or more \link cvm::atom_group +/// apply forces to atoms. Typically, one or more \link colvarmodule::atom_group /// \endlink objects are used, but this is not a requirement for as long as -/// the \link cvc \endlink object communicates with the simulation program. +/// the \link colvar::cvc \endlink object communicates with the simulation program. /// /// If you wish to implement a new collective variable component, you /// should write your own class by inheriting directly from \link @@ -75,9 +75,9 @@ public: /// \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 + /// object within its constructor, when all \link colvar::cvc \endlink /// objects are initialized; therefore the main "config string" - /// constructor does not need to define it. If a \link cvc + /// constructor does not need to define it. If a \link colvar::cvc /// \endlink is initialized and/or a different constructor is used, /// this variable definition should be set within the constructor. std::string function_type; @@ -109,6 +109,9 @@ public: /// cvc \endlink virtual int init(std::string const &conf); + /// \brief Initialize dependency tree + virtual int init_dependencies(); + /// \brief Within the constructor, make a group parse its own /// options from the provided configuration string /// Returns reference to new group @@ -122,7 +125,7 @@ public: /// \brief After construction, set data related to dependency handling int setup(); - /// \brief Default constructor (used when \link cvc \endlink + /// \brief Default constructor (used when \link colvar::cvc \endlink /// objects are declared within other ones) cvc(); @@ -133,7 +136,7 @@ public: static std::vector cvc_features; /// \brief Implementation of the feature list accessor for colvar - virtual const std::vector &features() + virtual const std::vector &features() const { return cvc_features; } @@ -148,6 +151,9 @@ public: cvc_features.clear(); } + /// \brief Get vector of vectors of atom IDs for all atom groups + virtual std::vector > get_atom_lists(); + /// \brief Obtain data needed for the calculation for the backend virtual void read_data(); @@ -164,6 +170,10 @@ public: /// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component virtual void debug_gradients(); + /// \brief Calculate atomic gradients and add them to the corresponding item in gradient vector + /// May be overridden by CVCs that do not store their gradients in the classic way, see dihedPC + virtual void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); + /// \brief Calculate the total force from the system using the /// inverse atomic gradients virtual void calc_force_invgrads(); @@ -237,7 +247,7 @@ public: colvarvalue const &x2) const; /// \brief Wrap value (for periodic/symmetric cvcs) - virtual void wrap(colvarvalue &x) const; + virtual void wrap(colvarvalue &x_unwrapped) const; /// \brief Pointers to all atom groups, to let colvars collect info /// e.g. atomic gradients @@ -246,7 +256,7 @@ public: /// \brief Store a pointer to new atom group, and list as child for dependencies inline void register_atom_group(cvm::atom_group *ag) { atom_groups.push_back(ag); - add_child((colvardeps *) ag); + add_child(ag); } /// \brief Whether or not this CVC will be computed in parallel whenever possible @@ -415,7 +425,7 @@ public: virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; /// \brief Redefined to make use of the user-provided period - virtual void wrap(colvarvalue &x) const; + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -474,7 +484,7 @@ public: virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x) const; + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -559,6 +569,35 @@ public: +/// \brief Colvar component: dipole magnitude of a molecule +class colvar::dipole_magnitude + : public colvar::cvc +{ +protected: + /// Dipole atom group + cvm::atom_group *atoms; + cvm::atom_pos dipoleV; +public: + /// Initialize by parsing the configuration + dipole_magnitude (std::string const &conf); + dipole_magnitude (cvm::atom const &a1); + dipole_magnitude(); + virtual inline ~dipole_magnitude() {} + 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; +}; + + + /// \brief Colvar component: Radius of gyration of an atom group /// (colvarvalue::type_scalar type, range [0:*)) class colvar::gyration @@ -818,7 +857,7 @@ public: virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x) const; + virtual void wrap(colvarvalue &x_unwrapped) const; }; @@ -1002,7 +1041,7 @@ public: cvm::atom const &donor, cvm::real r0, int en, int ed); h_bond(); - virtual ~h_bond(); + virtual ~h_bond() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); @@ -1090,6 +1129,8 @@ public: virtual ~alpha_angles(); void calc_value(); void calc_gradients(); + /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs + void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); void apply_force(colvarvalue const &force); virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; @@ -1120,6 +1161,8 @@ public: virtual ~dihedPC(); void calc_value(); void calc_gradients(); + /// Re-implementation of cvc::collect_gradients() to carry over atomic gradients of sub-cvcs + void collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients); void apply_force(colvarvalue const &force); virtual cvm::real dist2(colvarvalue const &x1, colvarvalue const &x2) const; @@ -1159,6 +1202,7 @@ public: orientation(std::string const &conf); orientation(); + virtual int init(std::string const &conf); virtual ~orientation() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1183,6 +1227,7 @@ public: orientation_angle(std::string const &conf); orientation_angle(); + virtual int init(std::string const &conf); virtual ~orientation_angle() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1207,6 +1252,7 @@ public: orientation_proj(std::string const &conf); orientation_proj(); + virtual int init(std::string const &conf); virtual ~orientation_proj() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1234,6 +1280,7 @@ public: tilt(std::string const &conf); tilt(); + virtual int init(std::string const &conf); virtual ~tilt() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1261,6 +1308,7 @@ public: spin_angle(std::string const &conf); spin_angle(); + virtual int init(std::string const &conf); virtual ~spin_angle() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1275,7 +1323,7 @@ public: virtual colvarvalue dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const; /// Redefined to handle the 2*PI periodicity - virtual void wrap(colvarvalue &x) const; + virtual void wrap(colvarvalue &x_unwrapped) const; }; diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 9f879a4c41..97fb23b181 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -11,9 +11,6 @@ #include "colvar.h" #include "colvarcomp.h" -#include - - colvar::angle::angle(std::string const &conf) : cvc(conf) @@ -77,14 +74,14 @@ void colvar::angle::calc_value() cvm::real const cos_theta = (r21*r23)/(r21l*r23l); - x.real_value = (180.0/PI) * std::acos(cos_theta); + x.real_value = (180.0/PI) * cvm::acos(cos_theta); } void colvar::angle::calc_gradients() { cvm::real const cos_theta = (r21*r23)/(r21l*r23l); - cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta); + cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta); dxdr1 = (180.0/PI) * dxdcos * (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l ); @@ -126,7 +123,7 @@ void colvar::angle::calc_Jacobian_derivative() // det(J) = (2 pi) r^2 * sin(theta) // hence Jd = cot(theta) const cvm::real theta = x.real_value * PI / 180.0; - jd = PI / 180.0 * (theta != 0.0 ? std::cos(theta) / std::sin(theta) : 0.0); + jd = PI / 180.0 * (theta != 0.0 ? cvm::cos(theta) / cvm::sin(theta) : 0.0); } @@ -202,7 +199,7 @@ void colvar::dipole_angle::calc_value() cvm::real const cos_theta = (r21*r23)/(r21l*r23l); - x.real_value = (180.0/PI) * std::acos(cos_theta); + x.real_value = (180.0/PI) * cvm::acos(cos_theta); } //to be implemented @@ -212,7 +209,7 @@ void colvar::dipole_angle::calc_value() void colvar::dipole_angle::calc_gradients() { cvm::real const cos_theta = (r21*r23)/(r21l*r23l); - cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta); + cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta); dxdr1 = (180.0/PI) * dxdcos * (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l ); @@ -346,7 +343,7 @@ void colvar::dihedral::calc_value() cvm::real const cos_phi = n1 * n2; cvm::real const sin_phi = n1 * r34 * r23.norm(); - x.real_value = (180.0/PI) * std::atan2(sin_phi, cos_phi); + x.real_value = (180.0/PI) * cvm::atan2(sin_phi, cos_phi); this->wrap(x); } @@ -368,7 +365,7 @@ void colvar::dihedral::calc_gradients() rB = 1.0/rB; B *= rB; - if (std::fabs(sin_phi) > 0.1) { + if (cvm::fabs(sin_phi) > 0.1) { rA = 1.0/rA; A *= rA; cvm::rvector const dcosdA = rA*(cos_phi*A-B); @@ -440,8 +437,8 @@ void colvar::dihedral::calc_force_invgrads() cvm::real const dot1 = u23 * u12; cvm::real const dot4 = u23 * u34; - cvm::real const fact1 = d12 * std::sqrt(1.0 - dot1 * dot1); - cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4); + cvm::real const fact1 = d12 * cvm::sqrt(1.0 - dot1 * dot1); + cvm::real const fact4 = d34 * cvm::sqrt(1.0 - dot4 * dot4); group1->read_total_forces(); if (is_enabled(f_cvc_one_site_total_force)) { @@ -508,19 +505,17 @@ colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1, } -void colvar::dihedral::wrap(colvarvalue &x) const +void colvar::dihedral::wrap(colvarvalue &x_unwrapped) const { - if ((x.real_value - wrap_center) >= 180.0) { - x.real_value -= 360.0; + if ((x_unwrapped.real_value - wrap_center) >= 180.0) { + x_unwrapped.real_value -= 360.0; return; } - if ((x.real_value - wrap_center) < -180.0) { - x.real_value += 360.0; + if ((x_unwrapped.real_value - wrap_center) < -180.0) { + x_unwrapped.real_value += 360.0; return; } - - return; } @@ -548,8 +543,8 @@ void colvar::polar_theta::calc_value() cvm::rvector pos = atoms->center_of_mass(); r = atoms->center_of_mass().norm(); // Internal values of theta and phi are radians - theta = (r > 0.) ? std::acos(pos.z / r) : 0.; - phi = std::atan2(pos.y, pos.x); + theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.; + phi = cvm::atan2(pos.y, pos.x); x.real_value = (180.0/PI) * theta; } @@ -560,9 +555,9 @@ void colvar::polar_theta::calc_gradients() atoms->set_weighted_gradient(cvm::rvector(0., 0., 0.)); else atoms->set_weighted_gradient(cvm::rvector( - (180.0/PI) * std::cos(theta) * std::cos(phi) / r, - (180.0/PI) * std::cos(theta) * std::sin(phi) / r, - (180.0/PI) * -std::sin(theta) / r)); + (180.0/PI) * cvm::cos(theta) * cvm::cos(phi) / r, + (180.0/PI) * cvm::cos(theta) * cvm::sin(phi) / r, + (180.0/PI) * -cvm::sin(theta) / r)); } @@ -602,8 +597,8 @@ void colvar::polar_phi::calc_value() cvm::rvector pos = atoms->center_of_mass(); r = atoms->center_of_mass().norm(); // Internal values of theta and phi are radians - theta = (r > 0.) ? std::acos(pos.z / r) : 0.; - phi = std::atan2(pos.y, pos.x); + theta = (r > 0.) ? cvm::acos(pos.z / r) : 0.; + phi = cvm::atan2(pos.y, pos.x); x.real_value = (180.0/PI) * phi; } @@ -611,8 +606,8 @@ void colvar::polar_phi::calc_value() void colvar::polar_phi::calc_gradients() { atoms->set_weighted_gradient(cvm::rvector( - (180.0/PI) * -std::sin(phi) / (r*std::sin(theta)), - (180.0/PI) * std::cos(phi) / (r*std::sin(theta)), + (180.0/PI) * -cvm::sin(phi) / (r*cvm::sin(theta)), + (180.0/PI) * cvm::cos(phi) / (r*cvm::sin(theta)), 0.)); } @@ -653,15 +648,15 @@ colvarvalue colvar::polar_phi::dist2_rgrad(colvarvalue const &x1, } -void colvar::polar_phi::wrap(colvarvalue &x) const +void colvar::polar_phi::wrap(colvarvalue &x_unwrapped) const { - if ((x.real_value - wrap_center) >= 180.0) { - x.real_value -= 360.0; + if ((x_unwrapped.real_value - wrap_center) >= 180.0) { + x_unwrapped.real_value -= 360.0; return; } - if ((x.real_value - wrap_center) < -180.0) { - x.real_value += 360.0; + if ((x_unwrapped.real_value - wrap_center) < -180.0) { + x_unwrapped.real_value += 360.0; return; } diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index ec53391ef5..059b0c825b 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -7,8 +7,6 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include - #include "colvarmodule.h" #include "colvarparse.h" #include "colvaratoms.h" @@ -102,6 +100,12 @@ colvar::coordnum::coordnum(std::string const &conf) group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); + if (group1 == NULL || group2 == NULL) { + cvm::error("Error: failed to initialize atom groups.\n", + INPUT_ERROR); + return; + } + if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) { cvm::error("Error: group1 and group2 share a common atom (number: " + cvm::to_str(atom_number) + ")\n", INPUT_ERROR); @@ -408,12 +412,6 @@ colvar::h_bond::h_bond() } -colvar::h_bond::~h_bond() -{ - delete atom_groups[0]; -} - - void colvar::h_bond::calc_value() { int const flags = coordnum::ef_null; @@ -655,8 +653,6 @@ colvar::groupcoordnum::groupcoordnum() void colvar::groupcoordnum::calc_value() { - cvm::rvector const r0_vec(0.0); // TODO enable the flag? - // create fake atoms to hold the com coordinates cvm::atom group1_com_atom; cvm::atom group2_com_atom; @@ -680,8 +676,6 @@ void colvar::groupcoordnum::calc_value() void colvar::groupcoordnum::calc_gradients() { - cvm::rvector const r0_vec(0.0); // TODO enable the flag? - cvm::atom group1_com_atom; cvm::atom group2_com_atom; group1_com_atom.pos = group1->center_of_mass(); diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index a2b1a5cd23..d9cd9d55e4 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -7,8 +7,6 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include - #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" @@ -102,7 +100,7 @@ colvar::distance_vec::distance_vec(std::string const &conf) { function_type = "distance_vec"; enable(f_cvc_com_based); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_3vector); } @@ -112,7 +110,7 @@ colvar::distance_vec::distance_vec() { function_type = "distance_vec"; enable(f_cvc_com_based); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_3vector); } @@ -320,7 +318,7 @@ cvm::real colvar::distance_z::dist2(colvarvalue const &x1, { cvm::real diff = x1.real_value - x2.real_value; if (b_periodic) { - cvm::real shift = std::floor(diff/period + 0.5); + cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } return diff * diff; @@ -332,7 +330,7 @@ colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1, { cvm::real diff = x1.real_value - x2.real_value; if (b_periodic) { - cvm::real shift = std::floor(diff/period + 0.5); + cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } return 2.0 * diff; @@ -344,22 +342,23 @@ colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1, { cvm::real diff = x1.real_value - x2.real_value; if (b_periodic) { - cvm::real shift = std::floor(diff/period + 0.5); + cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } return (-2.0) * diff; } -void colvar::distance_z::wrap(colvarvalue &x) const +void colvar::distance_z::wrap(colvarvalue &x_unwrapped) const { if (!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; + cvm::real shift = + cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5); + x_unwrapped.real_value -= shift * period; return; } @@ -481,7 +480,7 @@ colvar::distance_dir::distance_dir(std::string const &conf) { function_type = "distance_dir"; enable(f_cvc_com_based); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_unit3vector); } @@ -491,7 +490,7 @@ colvar::distance_dir::distance_dir() { function_type = "distance_dir"; enable(f_cvc_com_based); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_unit3vector); } @@ -629,7 +628,7 @@ void colvar::distance_inv::calc_value() } x.real_value *= 1.0 / cvm::real(group1->size() * group2->size()); - x.real_value = std::pow(x.real_value, -1.0/cvm::real(exponent)); + x.real_value = cvm::pow(x.real_value, -1.0/cvm::real(exponent)); cvm::real const dxdsum = (-1.0/(cvm::real(exponent))) * cvm::integer_power(x.real_value, exponent+1) / @@ -671,7 +670,7 @@ colvar::distance_pairs::distance_pairs(std::string const &conf) group2 = parse_group(conf, "group2"); x.type(colvarvalue::type_vector); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.vector1d_value.resize(group1->size() * group2->size()); } @@ -679,7 +678,7 @@ colvar::distance_pairs::distance_pairs(std::string const &conf) colvar::distance_pairs::distance_pairs() { function_type = "distance_pairs"; - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_vector); } @@ -747,6 +746,63 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force) +colvar::dipole_magnitude::dipole_magnitude(std::string const &conf) + : cvc(conf) +{ + function_type = "dipole_magnitude"; + atoms = parse_group(conf, "atoms"); + init_total_force_params(conf); + x.type(colvarvalue::type_scalar); +} + + +colvar::dipole_magnitude::dipole_magnitude(cvm::atom const &a1) +{ + atoms = new cvm::atom_group(std::vector(1, a1)); + register_atom_group(atoms); + x.type(colvarvalue::type_scalar); +} + + +colvar::dipole_magnitude::dipole_magnitude() +{ + function_type = "dipole_magnitude"; + x.type(colvarvalue::type_scalar); +} + + +void colvar::dipole_magnitude::calc_value() +{ + cvm::atom_pos const atomsCom = atoms->center_of_mass(); + atoms->calc_dipole(atomsCom); + dipoleV = atoms->dipole(); + x.real_value = dipoleV.norm(); +} + + +void colvar::dipole_magnitude::calc_gradients() +{ + cvm::real const aux1 = atoms->total_charge/atoms->total_mass; + cvm::atom_pos const dipVunit = dipoleV.unit(); + + for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) { + ai->grad = (ai->charge - aux1*ai->mass) * dipVunit; + } +} + + +void colvar::dipole_magnitude::apply_force(colvarvalue const &force) +{ + if (!atoms->noforce) { + atoms->apply_colvar_force(force.real_value); + } +} + + +simple_scalar_dist_functions(dipole_magnitude) + + + colvar::gyration::gyration(std::string const &conf) : cvc(conf) { @@ -782,7 +838,7 @@ void colvar::gyration::calc_value() 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())); + x.real_value = cvm::sqrt(x.real_value / cvm::real(atoms->size())); } @@ -1029,7 +1085,7 @@ void colvar::rmsd::calc_value() x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2(); } x.real_value /= cvm::real(atoms->size()); // MSD - x.real_value = std::sqrt(x.real_value); + x.real_value = cvm::sqrt(x.real_value); } @@ -1405,7 +1461,7 @@ void colvar::eigenvector::calc_Jacobian_derivative() } } - jd.real_value = sum * std::sqrt(eigenvec_invnorm2); + jd.real_value = sum * cvm::sqrt(eigenvec_invnorm2); } @@ -1436,7 +1492,7 @@ colvar::cartesian::cartesian(std::string const &conf) } x.type(colvarvalue::type_vector); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.vector1d_value.resize(atoms->size() * axes.size()); } diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index 91e47f13d9..b9f9c60cdb 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -7,8 +7,6 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include - #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" @@ -27,7 +25,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) cvm::log("Initializing alpha_angles object.\n"); function_type = "alpha_angles"; - enable(f_cvc_implicit_gradient); + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); std::string segment_id; @@ -118,7 +116,7 @@ colvar::alpha_angles::alpha_angles() : cvc() { function_type = "alpha_angles"; - enable(f_cvc_implicit_gradient); + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -133,6 +131,8 @@ colvar::alpha_angles::~alpha_angles() delete hb.back(); hb.pop_back(); } + // Our references to atom groups have become invalid now that children cvcs are deleted + atom_groups.clear(); } @@ -191,6 +191,58 @@ void colvar::alpha_angles::calc_gradients() } +void colvar::alpha_angles::collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients) +{ + cvm::real cvc_coeff = sup_coeff * cvm::real(sup_np) * cvm::integer_power(value().real_value, sup_np-1); + + if (theta.size()) { + cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real(theta.size()); + + for (size_t i = 0; i < theta.size(); i++) { + cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol; + cvm::real const f = ( (1.0 - (t*t)) / + (1.0 - (t*t*t*t)) ); + cvm::real const dfdt = + 1.0/(1.0 - (t*t*t*t)) * + ( (-2.0 * t) + (-1.0*f)*(-4.0 * (t*t*t)) ); + + // Coeficient of this CVC's gradient in the colvar gradient, times coefficient of this + // angle's gradient in the CVC's gradient + cvm::real const coeff = cvc_coeff * theta_norm * dfdt * (1.0/theta_tol); + + for (size_t j = 0; j < theta[i]->atom_groups.size(); j++) { + cvm::atom_group &ag = *(theta[i]->atom_groups[j]); + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; + } + } + } + } + + if (hb.size()) { + + cvm::real const hb_norm = hb_coeff / cvm::real(hb.size()); + + for (size_t i = 0; i < hb.size(); i++) { + // Coeficient of this CVC's gradient in the colvar gradient, times coefficient of this + // hbond's gradient in the CVC's gradient + cvm::real const coeff = cvc_coeff * 0.5 * hb_norm; + + for (size_t j = 0; j < hb[i]->atom_groups.size(); j++) { + cvm::atom_group &ag = *(hb[i]->atom_groups[j]); + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; + } + } + } + } +} + + void colvar::alpha_angles::apply_force(colvarvalue const &force) { @@ -242,7 +294,8 @@ colvar::dihedPC::dihedPC(std::string const &conf) cvm::log("Initializing dihedral PC object.\n"); function_type = "dihedPC"; - enable(f_cvc_implicit_gradient); + // Supported through references to atom groups of children cvcs + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); std::string segment_id; @@ -372,7 +425,8 @@ colvar::dihedPC::dihedPC() : cvc() { function_type = "dihedPC"; - enable(f_cvc_implicit_gradient); + // Supported through references to atom groups of children cvcs + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -383,6 +437,8 @@ colvar::dihedPC::~dihedPC() delete theta.back(); theta.pop_back(); } + // Our references to atom groups have become invalid now that children cvcs are deleted + atom_groups.clear(); } @@ -392,8 +448,8 @@ void colvar::dihedPC::calc_value() for (size_t i = 0; i < theta.size(); i++) { theta[i]->calc_value(); cvm::real const t = (PI / 180.) * theta[i]->value().real_value; - x.real_value += coeffs[2*i ] * std::cos(t) - + coeffs[2*i+1] * std::sin(t); + x.real_value += coeffs[2*i ] * cvm::cos(t) + + coeffs[2*i+1] * cvm::sin(t); } } @@ -406,12 +462,35 @@ void colvar::dihedPC::calc_gradients() } +void colvar::dihedPC::collect_gradients(std::vector const &atom_ids, std::vector &atomic_gradients) +{ + cvm::real cvc_coeff = sup_coeff * cvm::real(sup_np) * cvm::integer_power(value().real_value, sup_np-1); + for (size_t i = 0; i < theta.size(); i++) { + cvm::real const t = (PI / 180.) * theta[i]->value().real_value; + cvm::real const dcosdt = - (PI / 180.) * cvm::sin(t); + cvm::real const dsindt = (PI / 180.) * cvm::cos(t); + // Coeficient of this dihedPC's gradient in the colvar gradient, times coefficient of this + // dihedral's gradient in the dihedPC's gradient + cvm::real const coeff = cvc_coeff * (coeffs[2*i] * dcosdt + coeffs[2*i+1] * dsindt); + + for (size_t j = 0; j < theta[i]->atom_groups.size(); j++) { + cvm::atom_group &ag = *(theta[i]->atom_groups[j]); + for (size_t k = 0; k < ag.size(); k++) { + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; + } + } + } +} + + void colvar::dihedPC::apply_force(colvarvalue const &force) { for (size_t i = 0; i < theta.size(); i++) { cvm::real const t = (PI / 180.) * theta[i]->value().real_value; - cvm::real const dcosdt = - (PI / 180.) * std::sin(t); - cvm::real const dsindt = (PI / 180.) * std::cos(t); + cvm::real const dcosdt = - (PI / 180.) * cvm::sin(t); + cvm::real const dsindt = (PI / 180.) * cvm::cos(t); theta[i]->apply_force((coeffs[2*i ] * dcosdt + coeffs[2*i+1] * dsindt) * force); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 498ef7c2f5..3c8d0b4af6 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -7,8 +7,6 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include - #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" @@ -18,21 +16,27 @@ colvar::orientation::orientation(std::string const &conf) - : cvc(conf) + : cvc() { function_type = "orientation"; - atoms = parse_group(conf, "atoms"); - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_quaternion); + init(conf); +} + +int colvar::orientation::init(std::string const &conf) +{ + int error_code = cvc::init(conf); + + atoms = parse_group(conf, "atoms"); ref_pos.reserve(atoms->size()); 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::error("Error: reference positions do not " - "match the number of requested atoms.\n"); - return; + return cvm::error("Error: reference positions do not " + "match the number of requested atoms.\n", INPUT_ERROR); } } @@ -46,9 +50,8 @@ colvar::orientation::orientation(std::string const &conf) // use PDB flags if column is provided bool found = get_keyval(conf, "refPositionsColValue", file_col_value, 0.0); if (found && file_col_value==0.0) { - cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); - return; + return cvm::error("Error: refPositionsColValue, " + "if provided, must be non-zero.\n", INPUT_ERROR); } } @@ -59,9 +62,8 @@ colvar::orientation::orientation(std::string const &conf) } if (!ref_pos.size()) { - cvm::error("Error: must define a set of " - "reference coordinates.\n"); - return; + return cvm::error("Error: must define a set of " + "reference coordinates.\n", INPUT_ERROR); } @@ -85,6 +87,7 @@ colvar::orientation::orientation(std::string const &conf) rot.request_group2_gradients(atoms->size()); } + return error_code; } @@ -92,7 +95,7 @@ colvar::orientation::orientation() : cvc() { function_type = "orientation"; - enable(f_cvc_implicit_gradient); + disable(f_cvc_explicit_gradient); x.type(colvarvalue::type_quaternion); } @@ -158,10 +161,18 @@ colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1, colvar::orientation_angle::orientation_angle(std::string const &conf) - : orientation(conf) + : orientation() { function_type = "orientation_angle"; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); + init(conf); +} + + +int colvar::orientation_angle::init(std::string const &conf) +{ + return orientation::init(conf); } @@ -169,6 +180,7 @@ colvar::orientation_angle::orientation_angle() : orientation() { function_type = "orientation_angle"; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -180,9 +192,9 @@ void colvar::orientation_angle::calc_value() rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); if ((rot.q).q0 >= 0.0) { - x.real_value = (180.0/PI) * 2.0 * std::acos((rot.q).q0); + x.real_value = (180.0/PI) * 2.0 * cvm::acos((rot.q).q0); } else { - x.real_value = (180.0/PI) * 2.0 * std::acos(-1.0 * (rot.q).q0); + x.real_value = (180.0/PI) * 2.0 * cvm::acos(-1.0 * (rot.q).q0); } } @@ -191,7 +203,7 @@ void colvar::orientation_angle::calc_gradients() { cvm::real const dxdq0 = ( ((rot.q).q0 * (rot.q).q0 < 1.0) ? - ((180.0 / PI) * (-2.0) / std::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) : + ((180.0 / PI) * (-2.0) / cvm::sqrt(1.0 - ((rot.q).q0 * (rot.q).q0))) : 0.0 ); for (size_t ia = 0; ia < atoms->size(); ia++) { @@ -214,10 +226,18 @@ simple_scalar_dist_functions(orientation_angle) colvar::orientation_proj::orientation_proj(std::string const &conf) - : orientation(conf) + : orientation() { function_type = "orientation_proj"; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); + init(conf); +} + + +int colvar::orientation_proj::init(std::string const &conf) +{ + return orientation::init(conf); } @@ -225,6 +245,7 @@ colvar::orientation_proj::orientation_proj() : orientation() { function_type = "orientation_proj"; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -261,18 +282,28 @@ simple_scalar_dist_functions(orientation_proj) colvar::tilt::tilt(std::string const &conf) - : orientation(conf) + : orientation() { function_type = "tilt"; + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); + init(conf); +} + + +int colvar::tilt::init(std::string const &conf) +{ + int error_code = COLVARS_OK; + + error_code |= orientation::init(conf); get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0)); - if (axis.norm2() != 1.0) { axis /= axis.norm(); cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n"); } - x.type(colvarvalue::type_scalar); + return error_code; } @@ -280,6 +311,7 @@ colvar::tilt::tilt() : orientation() { function_type = "tilt"; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -322,20 +354,30 @@ simple_scalar_dist_functions(tilt) colvar::spin_angle::spin_angle(std::string const &conf) - : orientation(conf) + : orientation() { function_type = "spin_angle"; + period = 360.0; + b_periodic = true; + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); + init(conf); +} + + +int colvar::spin_angle::init(std::string const &conf) +{ + int error_code = COLVARS_OK; + + error_code |= orientation::init(conf); get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0)); - if (axis.norm2() != 1.0) { axis /= axis.norm(); cvm::log("Normalizing rotation axis to "+cvm::to_str(axis)+".\n"); } - period = 360.0; - b_periodic = true; - x.type(colvarvalue::type_scalar); + return error_code; } @@ -345,6 +387,7 @@ colvar::spin_angle::spin_angle() function_type = "spin_angle"; period = 360.0; b_periodic = true; + enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } @@ -410,15 +453,15 @@ colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1, } -void colvar::spin_angle::wrap(colvarvalue &x) const +void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const { - if ((x.real_value - wrap_center) >= 180.0) { - x.real_value -= 360.0; + if ((x_unwrapped.real_value - wrap_center) >= 180.0) { + x_unwrapped.real_value -= 360.0; return; } - if ((x.real_value - wrap_center) < -180.0) { - x.real_value += 360.0; + if ((x_unwrapped.real_value - wrap_center) < -180.0) { + x_unwrapped.real_value += 360.0; return; } diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index d20ee6e55c..276f2b39e7 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -122,12 +122,7 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp, int colvardeps::enable(int feature_id, bool dry_run /* default: false */, - // dry_run: fail silently, do not enable if available - // flag is passed recursively to deps of this feature bool toplevel /* default: true */) -// toplevel: false if this is called as part of a chain of dependency resolution -// this is used to diagnose failed dependencies by displaying the full stack -// only the toplevel dependency will throw a fatal error { int res; size_t i, j; @@ -143,8 +138,7 @@ int colvardeps::enable(int feature_id, if (fs->enabled) { if (!(dry_run || toplevel)) { - // This is a dependency - // Prevent disabling this feature as long + // This is a dependency: prevent disabling this feature as long // as requirement is enabled fs->ref_count++; if (cvm::debug()) @@ -173,7 +167,10 @@ int colvardeps::enable(int feature_id, if (!toplevel && !is_dynamic(feature_id)) { if (!dry_run) { cvm::log(feature_type_descr + " feature \"" + f->description - + "\" may not be enabled as a dependency in " + description + ".\n"); + + "\" cannot be enabled automatically in " + description + "."); + if (is_user(feature_id)) { + cvm::log("Try setting it manually.\n"); + } } return COLVARS_ERROR; } @@ -354,6 +351,7 @@ int colvardeps::disable(int feature_id) { return COLVARS_OK; } + int colvardeps::decr_ref_count(int feature_id) { int &rc = feature_states[feature_id].ref_count; feature *f = features()[feature_id]; @@ -379,324 +377,52 @@ int colvardeps::decr_ref_count(int feature_id) { return COLVARS_OK; } -void colvardeps::init_feature(int feature_id, const char *description, feature_type type) { - modify_features()[feature_id]->description = description; + +void colvardeps::init_feature(int feature_id, const char *description_in, feature_type type) { + modify_features()[feature_id]->description = description_in; modify_features()[feature_id]->type = type; } -// Shorthand macros for describing dependencies -#define f_req_self(f, g) features()[f]->requires_self.push_back(g) -// This macro ensures that exclusions are symmetric -#define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \ - features()[g]->requires_exclude.push_back(f) -#define f_req_children(f, g) features()[f]->requires_children.push_back(g) -#define f_req_alt2(f, g, h) features()[f]->requires_alt.push_back(std::vector(2));\ - features()[f]->requires_alt.back()[0] = g; \ - features()[f]->requires_alt.back()[1] = h -#define f_req_alt3(f, g, h, i) features()[f]->requires_alt.push_back(std::vector(3));\ - features()[f]->requires_alt.back()[0] = g; \ - features()[f]->requires_alt.back()[1] = h; \ - features()[f]->requires_alt.back()[2] = i -#define f_req_alt4(f, g, h, i, j) features()[f]->requires_alt.push_back(std::vector(4));\ - features()[f]->requires_alt.back()[0] = g; \ - features()[f]->requires_alt.back()[1] = h; \ - features()[f]->requires_alt.back()[2] = i; \ - features()[f]->requires_alt.back()[3] = j -void colvardeps::init_cvb_requires() { - int i; - if (features().size() == 0) { - for (i = 0; i < f_cvb_ntot; i++) { - modify_features().push_back(new feature); - } - - init_feature(f_cvb_active, "active", f_type_dynamic); - f_req_children(f_cvb_active, f_cv_active); - - init_feature(f_cvb_awake, "awake", f_type_static); - f_req_self(f_cvb_awake, f_cvb_active); - - init_feature(f_cvb_apply_force, "apply force", f_type_user); - f_req_children(f_cvb_apply_force, f_cv_gradient); - - init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic); - f_req_children(f_cvb_get_total_force, f_cv_total_force); - - init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user); - f_req_self(f_cvb_output_acc_work, f_cvb_apply_force); - - init_feature(f_cvb_history_dependent, "history-dependent", f_type_static); - - init_feature(f_cvb_time_dependent, "time-dependent", f_type_static); - - init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static); - f_req_children(f_cvb_scalar_variables, f_cv_scalar); - - init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static); - - init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic); - f_req_self(f_cvb_calc_ti_samples, f_cvb_get_total_force); - f_req_children(f_cvb_calc_ti_samples, f_cv_grid); - - init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user); - f_req_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples); - - init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user); - f_req_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples); - } - - // Initialize feature_states for each instance - feature_states.reserve(f_cvb_ntot); - for (i = 0; i < f_cvb_ntot; i++) { - feature_states.push_back(feature_state(true, false)); - // Most features are available, so we set them so - // and list exceptions below - } - - // only compute TI samples when deriving from colvarbias_ti - feature_states[f_cvb_calc_ti_samples].available = false; +// Shorthand functions for describing dependencies +void colvardeps::require_feature_self(int f, int g) { + features()[f]->requires_self.push_back(g); } -void colvardeps::init_cv_requires() { - size_t i; - if (features().size() == 0) { - for (i = 0; i < f_cv_ntot; i++) { - modify_features().push_back(new feature); - } - - init_feature(f_cv_active, "active", f_type_dynamic); - // Do not require f_cvc_active in children, as some components may be disabled - // Colvars must be either a linear combination, or scalar (and polynomial) or scripted/custom - f_req_alt4(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted, f_cv_custom_function); - - init_feature(f_cv_awake, "awake", f_type_static); - f_req_self(f_cv_awake, f_cv_active); - - init_feature(f_cv_gradient, "gradient", f_type_dynamic); - f_req_children(f_cv_gradient, f_cvc_gradient); - - init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic); - f_req_self(f_cv_collect_gradient, f_cv_gradient); - f_req_self(f_cv_collect_gradient, f_cv_scalar); - // The following exlusion could be lifted by implementing the feature - f_req_exclude(f_cv_collect_gradient, f_cv_scripted); - - init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic); - - // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly - init_feature(f_cv_total_force, "total force", f_type_dynamic); - f_req_alt2(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc); - - // Deps for explicit total force calculation - init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic); - f_req_self(f_cv_total_force_calc, f_cv_scalar); - f_req_self(f_cv_total_force_calc, f_cv_linear); - f_req_children(f_cv_total_force_calc, f_cvc_inv_gradient); - f_req_self(f_cv_total_force_calc, f_cv_Jacobian); - - init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic); - f_req_self(f_cv_Jacobian, f_cv_scalar); - f_req_self(f_cv_Jacobian, f_cv_linear); - f_req_children(f_cv_Jacobian, f_cvc_Jacobian); - - init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user); - f_req_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated - - init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user); - f_req_self(f_cv_extended_Lagrangian, f_cv_scalar); - f_req_self(f_cv_extended_Lagrangian, f_cv_gradient); - - init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user); - f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian); - - init_feature(f_cv_linear, "linear", f_type_static); - - init_feature(f_cv_scalar, "scalar", f_type_static); - - init_feature(f_cv_output_energy, "output energy", f_type_user); - - init_feature(f_cv_output_value, "output value", f_type_user); - - init_feature(f_cv_output_velocity, "output velocity", f_type_user); - f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity); - - init_feature(f_cv_output_applied_force, "output applied force", f_type_user); - - init_feature(f_cv_output_total_force, "output total force", f_type_user); - f_req_self(f_cv_output_total_force, f_cv_total_force); - - init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user); - f_req_self(f_cv_subtract_applied_force, f_cv_total_force); - - init_feature(f_cv_lower_boundary, "lower boundary", f_type_user); - f_req_self(f_cv_lower_boundary, f_cv_scalar); - - init_feature(f_cv_upper_boundary, "upper boundary", f_type_user); - f_req_self(f_cv_upper_boundary, f_cv_scalar); - - init_feature(f_cv_grid, "grid", f_type_dynamic); - f_req_self(f_cv_grid, f_cv_lower_boundary); - f_req_self(f_cv_grid, f_cv_upper_boundary); - - init_feature(f_cv_runave, "running average", f_type_user); - - init_feature(f_cv_corrfunc, "correlation function", f_type_user); - - init_feature(f_cv_scripted, "scripted", f_type_user); - - init_feature(f_cv_custom_function, "custom function", f_type_user); - f_req_exclude(f_cv_custom_function, f_cv_scripted); - - init_feature(f_cv_periodic, "periodic", f_type_static); - f_req_self(f_cv_periodic, f_cv_scalar); - init_feature(f_cv_scalar, "scalar", f_type_static); - init_feature(f_cv_linear, "linear", f_type_static); - init_feature(f_cv_homogeneous, "homogeneous", f_type_static); - - // because total forces are obtained from the previous time step, - // we cannot (currently) have colvar values and total forces for the same timestep - init_feature(f_cv_multiple_ts, "multiple timestep colvar"); - f_req_exclude(f_cv_multiple_ts, f_cv_total_force_calc); - } - - // Initialize feature_states for each instance - feature_states.reserve(f_cv_ntot); - for (i = 0; i < f_cv_ntot; i++) { - feature_states.push_back(feature_state(true, false)); - // Most features are available, so we set them so - // and list exceptions below - } - - feature_states[f_cv_fdiff_velocity].available = - cvm::main()->proxy->simulation_running(); +// Ensure that exclusions are symmetric +void colvardeps::exclude_feature_self(int f, int g) { + features()[f]->requires_exclude.push_back(g); + features()[g]->requires_exclude.push_back(f); } -void colvardeps::init_cvc_requires() { - size_t i; - // Initialize static array once and for all - if (features().size() == 0) { - for (i = 0; i < colvardeps::f_cvc_ntot; i++) { - modify_features().push_back(new feature); - } - - init_feature(f_cvc_active, "active", f_type_dynamic); -// The dependency below may become useful if we use dynamic atom groups -// f_req_children(f_cvc_active, f_ag_active); - - init_feature(f_cvc_scalar, "scalar", f_type_static); - - init_feature(f_cvc_gradient, "gradient", f_type_dynamic); - - init_feature(f_cvc_implicit_gradient, "implicit gradient", f_type_static); - f_req_children(f_cvc_implicit_gradient, f_ag_implicit_gradient); - - init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic); - f_req_self(f_cvc_inv_gradient, f_cvc_gradient); - - init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user); - f_req_self(f_cvc_debug_gradient, f_cvc_gradient); - f_req_exclude(f_cvc_debug_gradient, f_cvc_implicit_gradient); - - init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic); - f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient); - - init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static); - - // init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user); - - // Compute total force on first site only to avoid unwanted - // coupling to other colvars (see e.g. Ciccotti et al., 2005) - init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user); - f_req_self(f_cvc_one_site_total_force, f_cvc_com_based); - - init_feature(f_cvc_scalable, "scalable calculation", f_type_static); - f_req_self(f_cvc_scalable, f_cvc_scalable_com); - - init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static); - f_req_self(f_cvc_scalable_com, f_cvc_com_based); - - - // TODO only enable this when f_ag_scalable can be turned on for a pre-initialized group - // f_req_children(f_cvc_scalable, f_ag_scalable); - // f_req_children(f_cvc_scalable_com, f_ag_scalable_com); - } - - // Initialize feature_states for each instance - // default as available, not enabled - // except dynamic features which default as unavailable - feature_states.reserve(f_cvc_ntot); - for (i = 0; i < colvardeps::f_cvc_ntot; i++) { - bool avail = is_dynamic(i) ? false : true; - feature_states.push_back(feature_state(avail, false)); - } - - // CVCs are enabled from the start - get disabled based on flags - feature_states[f_cvc_active].enabled = true; - - // Features that are implemented by all cvcs by default - // Each cvc specifies what other features are available - feature_states[f_cvc_active].available = true; - feature_states[f_cvc_gradient].available = true; - - // Use minimum-image distances by default - feature_states[f_cvc_pbc_minimum_image].enabled = true; - - // Features that are implemented by default if their requirements are - feature_states[f_cvc_one_site_total_force].available = true; - - // Features That are implemented only for certain simulation engine configurations - feature_states[f_cvc_scalable_com].available = (cvm::proxy->scalable_group_coms() == COLVARS_OK); - feature_states[f_cvc_scalable].available = feature_states[f_cvc_scalable_com].available; +void colvardeps::require_feature_children(int f, int g) { + features()[f]->requires_children.push_back(g); } -void colvardeps::init_ag_requires() { - size_t i; - // Initialize static array once and for all - if (features().size() == 0) { - for (i = 0; i < f_ag_ntot; i++) { - modify_features().push_back(new feature); - } +void colvardeps::require_feature_alt(int f, int g, int h) { + features()[f]->requires_alt.push_back(std::vector(2)); + features()[f]->requires_alt.back()[0] = g; + features()[f]->requires_alt.back()[1] = h; +} - init_feature(f_ag_active, "active", f_type_dynamic); - init_feature(f_ag_center, "translational fit", f_type_static); - init_feature(f_ag_rotate, "rotational fit", f_type_static); - init_feature(f_ag_fitting_group, "reference positions group", f_type_static); - init_feature(f_ag_implicit_gradient, "implicit atom gradient", f_type_dynamic); - init_feature(f_ag_fit_gradients, "fit gradients", f_type_user); - f_req_exclude(f_ag_fit_gradients, f_ag_implicit_gradient); - init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic); +void colvardeps::require_feature_alt(int f, int g, int h, int i) { + features()[f]->requires_alt.push_back(std::vector(3)); + features()[f]->requires_alt.back()[0] = g; + features()[f]->requires_alt.back()[1] = h; + features()[f]->requires_alt.back()[2] = i; +} - // parallel calculation implies that we have at least a scalable center of mass, - // but f_ag_scalable is kept as a separate feature to deal with future dependencies - init_feature(f_ag_scalable, "scalable group calculation", f_type_static); - init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static); - f_req_self(f_ag_scalable, f_ag_scalable_com); -// init_feature(f_ag_min_msd_fit, "minimum MSD fit") -// f_req_self(f_ag_min_msd_fit, f_ag_center) -// f_req_self(f_ag_min_msd_fit, f_ag_rotate) -// f_req_exclude(f_ag_min_msd_fit, f_ag_fitting_group) - } - - // Initialize feature_states for each instance - // default as unavailable, not enabled - feature_states.reserve(f_ag_ntot); - for (i = 0; i < colvardeps::f_ag_ntot; i++) { - feature_states.push_back(feature_state(false, false)); - } - - // Features that are implemented (or not) by all atom groups - feature_states[f_ag_active].available = true; - // f_ag_scalable_com is provided by the CVC iff it is COM-based - feature_states[f_ag_scalable_com].available = false; - // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else) - feature_states[f_ag_scalable].available = true; - feature_states[f_ag_fit_gradients].available = true; - feature_states[f_ag_implicit_gradient].available = true; +void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) { + features()[f]->requires_alt.push_back(std::vector(4)); + features()[f]->requires_alt.back()[0] = g; + features()[f]->requires_alt.back()[1] = h; + features()[f]->requires_alt.back()[2] = i; + features()[f]->requires_alt.back()[3] = j; } @@ -720,7 +446,7 @@ void colvardeps::print_state() { void colvardeps::add_child(colvardeps *child) { children.push_back(child); - child->parents.push_back((colvardeps *)this); + child->parents.push_back(this); // Solve dependencies of already enabled parent features // in the new child diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h index 940eefb01b..198a24f330 100644 --- a/lib/colvars/colvardeps.h +++ b/lib/colvars/colvardeps.h @@ -23,7 +23,11 @@ /// 3. Static features are static properties of the object, determined /// programatically at initialization time. /// -/// In all classes, feature 0 is active. When an object is inactivated +/// The following diagram summarizes the dependency tree at the bias, colvar, and colvarcomp levels. +/// Isolated and atom group features are not shown to save space. +/// @image html deps_2019.svg +/// +/// In all classes, feature 0 is `active`. When an object is inactivated /// all its children dependencies are dereferenced (free_children_deps) /// While the object is inactive, no dependency solving is done on children /// it is done when the object is activated back (restore_children_deps) @@ -72,7 +76,6 @@ protected: /// Unused by lower-level objects (cvcs and atom groups) int time_step_factor; -private: /// List of the states of all features std::vector feature_states; @@ -89,14 +92,14 @@ public: inline int get_time_step_factor() const {return time_step_factor;} /// Pair a numerical feature ID with a description and type - void init_feature(int feature_id, const char *description, feature_type type = f_type_not_set); + void init_feature(int feature_id, const char *description, feature_type type); /// Describes a feature and its dependencies /// used in a static array within each subclass class feature { public: - feature() {} + feature() : type(f_type_not_set) {} ~feature() {} std::string description; // Set by derived object initializer @@ -126,6 +129,7 @@ public: feature_type type; }; + inline bool is_not_set(int id) { return features()[id]->type == f_type_not_set; } inline bool is_dynamic(int id) { return features()[id]->type == f_type_dynamic; } inline bool is_static(int id) { return features()[id]->type == f_type_static; } inline bool is_user(int id) { return features()[id]->type == f_type_user; } @@ -135,7 +139,7 @@ public: // with a non-static array // Intermediate classes (colvarbias and colvarcomp, which are also base classes) // implement this as virtual to allow overriding - virtual const std::vector&features() = 0; + virtual const std::vector &features() const = 0; virtual std::vector&modify_features() = 0; void add_child(colvardeps *child); @@ -188,10 +192,12 @@ protected: public: - /// enable a feature and recursively solve its dependencies - /// for proper reference counting, one should not add - /// spurious calls to enable() - /// dry_run is set to true to recursively test if a feature is available, without enabling it + /// Enable a feature and recursively solve its dependencies. + /// For accurate reference counting, do not add spurious calls to enable() + /// \param dry_run Recursively test if a feature is available, without enabling it + /// \param toplevel False if this is called as part of a chain of dependency resolution. + /// This is used to diagnose failed dependencies by displaying the full stack: + /// only the toplevel dependency will throw a fatal error. int enable(int f, bool dry_run = false, bool toplevel = true); /// Disable a feature, decrease the reference count of its dependencies @@ -318,8 +324,8 @@ public: f_cvc_active, f_cvc_scalar, f_cvc_gradient, - /// \brief CVC doesn't calculate and store explicit atom gradients - f_cvc_implicit_gradient, + /// \brief CVC calculates and stores explicit atom gradients + f_cvc_explicit_gradient, f_cvc_inv_gradient, /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed f_cvc_debug_gradient, @@ -341,7 +347,7 @@ public: /// ie. not using refpositionsgroup // f_ag_min_msd_fit, /// \brief Does not have explicit atom gradients from parent CVC - f_ag_implicit_gradient, + f_ag_explicit_gradient, f_ag_fit_gradients, f_ag_atom_forces, f_ag_scalable, @@ -349,13 +355,39 @@ public: f_ag_ntot }; - void init_cvb_requires(); - void init_cv_requires(); - void init_cvc_requires(); - void init_ag_requires(); + /// Initialize dependency tree for object of a derived class + virtual int init_dependencies() = 0; + + /// Make feature f require feature g within the same object + void require_feature_self(int f, int g); + + /// Make features f and g mutually exclusive within the same object + void exclude_feature_self(int f, int g); + + /// Make feature f require feature g within children + void require_feature_children(int f, int g); + + /// Make feature f require either g or h within the same object + void require_feature_alt(int f, int g, int h); + + /// Make feature f require any of g, h, or i within the same object + void require_feature_alt(int f, int g, int h, int i); + + /// Make feature f require any of g, h, i, or j within the same object + void require_feature_alt(int f, int g, int h, int i, int j); /// \brief print all enabled features and those of children, for debugging void print_state(); + + /// \brief Check that a feature is enabled, raising BUG_ERROR if not + inline void check_enabled(int f, std::string const &reason) const + { + if (! is_enabled(f)) { + cvm::error("Error: "+reason+" requires that the feature \""+ + features()[f]->description+"\" is active.\n", BUG_ERROR); + } + } + }; #endif diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index 407b336afd..dc1a709edb 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -109,7 +109,7 @@ cvm::real colvar_grid_scalar::entropy() const { cvm::real sum = 0.0; for (size_t i = 0; i < nt; i++) { - sum += -1.0 * data[i] * std::log(data[i]); + sum += -1.0 * data[i] * cvm::logn(data[i]); } cvm::real bin_volume = 1.0; for (size_t id = 0; id < widths.size(); id++) { diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 9a0fe4c8ec..2ba0566e49 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -12,7 +12,6 @@ #include #include -#include #include "colvar.h" #include "colvarmodule.h" @@ -53,7 +52,7 @@ protected: std::vector cv; /// Do we request actual value (for extended-system colvars)? - std::vector actual_value; + std::vector use_actual_value; /// Get the low-level index corresponding to an index inline size_t address(std::vector const &ix) const @@ -136,8 +135,8 @@ public: inline void request_actual_value(bool b = true) { size_t i; - for (i = 0; i < actual_value.size(); i++) { - actual_value[i] = b; + for (i = 0; i < use_actual_value.size(); i++) { + use_actual_value[i] = b; } } @@ -215,7 +214,7 @@ public: mult(g.mult), data(), cv(g.cv), - actual_value(g.actual_value), + use_actual_value(g.use_actual_value), lower_boundaries(g.lower_boundaries), upper_boundaries(g.upper_boundaries), periodic(g.periodic), @@ -290,13 +289,13 @@ public: periodic.push_back(cv[i]->periodic_boundaries()); // By default, get reported colvar value (for extended Lagrangian colvars) - actual_value.push_back(false); + use_actual_value.push_back(false); // except if a colvar is specified twice in a row // then the first instance is the actual value // For histograms of extended-system coordinates if (i > 0 && cv[i-1] == cv[i]) { - actual_value[i-1] = true; + use_actual_value[i-1] = true; } if (margin) { @@ -319,8 +318,7 @@ public: return this->setup(); } - int init_from_boundaries(T const &t = T(), - size_t const &mult_i = 1) + int init_from_boundaries() { if (cvm::debug()) { cvm::log("Configuring grid dimensions from colvars boundaries.\n"); @@ -337,7 +335,7 @@ public: 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) { + if (cvm::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)+ @@ -392,20 +390,20 @@ public: /// \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); + return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); } /// \brief Report the bin corresponding to the current value of variable i /// and assign first or last bin if out of boundaries inline int current_bin_scalar_bound(int const i) const { - return value_to_bin_scalar_bound(actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); + return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); } /// \brief Report the bin corresponding to the current value of item iv in variable i inline int current_bin_scalar(int const i, int const iv) const { - return value_to_bin_scalar(actual_value[i] ? + return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value().vector1d_value[iv] : cv[i]->value().vector1d_value[iv], i); } @@ -414,14 +412,14 @@ public: /// 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] ); + return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); } /// \brief Use the lower boundary and the width to report which bin /// the provided value is in and assign first or last bin if out of boundaries inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const { - int bin_index = std::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); + int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); if (bin_index < 0) bin_index=0; if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1; return (int) bin_index; @@ -432,7 +430,7 @@ public: colvarvalue const &new_offset, cvm::real const &new_width) const { - return (int) std::floor( (value.real_value - new_offset.real_value) / new_width ); + return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width ); } /// \brief Use the two boundaries and the width to report the @@ -611,8 +609,8 @@ public: 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]; + cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i]; + cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i]; if (values[i].real_value < lower_boundaries[i]) dl *= -1.0; @@ -841,7 +839,7 @@ public: if (nd < lower_boundaries.size()) nd = lower_boundaries.size(); - if (! actual_value.size()) actual_value.assign(nd, false); + if (! use_actual_value.size()) use_actual_value.assign(nd, false); if (! periodic.size()) periodic.assign(nd, false); if (! widths.size()) widths.assign(nd, 1.0); @@ -849,7 +847,7 @@ public: if (old_nx.size()) { for (size_t i = 0; i < nd; i++) { if ( (old_nx[i] != nx[i]) || - (std::sqrt(cv[i]->dist2(old_lb[i], + (cvm::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > 1.0E-10) ) { new_params = true; } @@ -874,11 +872,11 @@ public: void check_consistency() { for (size_t i = 0; i < nd; i++) { - if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary, + if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary, lower_boundaries[i])) > 1.0E-10) || - (std::sqrt(cv[i]->dist2(cv[i]->upper_boundary, + (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary, upper_boundaries[i])) > 1.0E-10) || - (std::sqrt(cv[i]->dist2(cv[i]->width, + (cvm::sqrt(cv[i]->dist2(cv[i]->width, widths[i])) > 1.0E-10) ) { cvm::error("Error: restart information for a grid is " "inconsistent with that of its colvars.\n"); @@ -896,11 +894,11 @@ public: // 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] - + if ( (cvm::fabs(other_grid.lower_boundaries[i] - lower_boundaries[i]) > 1.0E-10) || - (std::fabs(other_grid.upper_boundaries[i] - + (cvm::fabs(other_grid.upper_boundaries[i] - upper_boundaries[i]) > 1.0E-10) || - (std::fabs(other_grid.widths[i] - + (cvm::fabs(other_grid.widths[i] - widths[i]) > 1.0E-10) || (data.size() != other_grid.data.size()) ) { cvm::error("Error: inconsistency between " @@ -1036,11 +1034,11 @@ public: std::istream & read_multicol(std::istream &is, bool add = false) { // Data in the header: nColvars, then for each - // xiMin, dXi, nPoints, periodic + // xiMin, dXi, nPoints, periodic flag std::string hash; cvm::real lower, width, x; - size_t n, periodic; + size_t n, periodic_flag; bool remap; std::vector new_value; std::vector nx_read; @@ -1053,7 +1051,8 @@ public: if ( !(is >> hash) || (hash != "#") ) { cvm::error("Error reading grid at position "+ - cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n"); + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); return is; } @@ -1075,15 +1074,16 @@ public: for (size_t i = 0; i < nd; i++ ) { if ( !(is >> hash) || (hash != "#") ) { cvm::error("Error reading grid at position "+ - cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n"); + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); return is; } - is >> lower >> width >> nx_read[i] >> periodic; + is >> lower >> width >> nx_read[i] >> periodic_flag; - if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) || - (std::fabs(width - widths[i] ) > 1.0e-10) || + if ( (cvm::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) || + (cvm::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"); @@ -1246,7 +1246,7 @@ public: if (A0 * A1 == 0) { return 0.; // can't handle empty bins } else { - return (std::log((cvm::real)A1) - std::log((cvm::real)A0)) + return (cvm::logn((cvm::real)A1) - cvm::logn((cvm::real)A0)) / (widths[n] * 2.); } } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge @@ -1258,7 +1258,7 @@ public: if (A0 * A1 == 0) { return 0.; // can't handle empty bins } else { - return (std::log((cvm::real)A1) - std::log((cvm::real)A0)) + return (cvm::logn((cvm::real)A1) - cvm::logn((cvm::real)A0)) / (widths[n] * 2.); } } else { @@ -1271,8 +1271,8 @@ public: if (A0 * A1 * A2 == 0) { return 0.; // can't handle empty bins } else { - return (-1.5 * std::log((cvm::real)A0) + 2. * std::log((cvm::real)A1) - - 0.5 * std::log((cvm::real)A2)) * increment / widths[n]; + return (-1.5 * cvm::logn((cvm::real)A0) + 2. * cvm::logn((cvm::real)A1) + - 0.5 * cvm::logn((cvm::real)A2)) * increment / widths[n]; } } } diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index d88a97a441..baffc25c28 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -28,6 +28,8 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) { depth_s = 0; + log_level_ = 10; + cv_traj_os = NULL; if (proxy == NULL) { @@ -152,12 +154,12 @@ int colvarmodule::read_config_string(std::string const &config_str) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration:\n"); - std::istringstream config_s(config_str); + std::istringstream new_config_s(config_str); // strip the comments away std::string conf = ""; std::string line; - while (parse->read_config_line(config_s, line)) { + while (parse->read_config_line(new_config_s, line)) { // Delete lines that contain only white space after removing comments if (line.find_first_not_of(colvarparse::white_space) != std::string::npos) conf.append(line+"\n"); @@ -255,13 +257,16 @@ int colvarmodule::append_new_config(std::string const &new_conf) int colvarmodule::parse_global_params(std::string const &conf) { - colvarmodule *cvm = cvm::main(); + // TODO document and then echo this keyword + parse->get_keyval(conf, "logLevel", log_level_, log_level_, + colvarparse::parse_silent); { std::string index_file_name; size_t pos = 0; while (parse->key_lookup(conf, "indexFile", &index_file_name, &pos)) { - cvm->read_index_file(index_file_name.c_str()); + cvm::log("# indexFile = \""+index_file_name+"\"\n"); + read_index_file(index_file_name.c_str()); index_file_name.clear(); } } @@ -580,6 +585,24 @@ cvm::atom_group *colvarmodule::atom_group_by_name(std::string const &name) } +void colvarmodule::register_named_atom_group(atom_group *ag) { + named_atom_groups.push_back(ag); +} + + +void colvarmodule::unregister_named_atom_group(cvm::atom_group *ag) +{ + for (std::vector::iterator agi = named_atom_groups.begin(); + agi != named_atom_groups.end(); + agi++) { + if (*agi == ag) { + named_atom_groups.erase(agi); + break; + } + } +} + + int colvarmodule::change_configuration(std::string const &bias_name, std::string const &conf) { @@ -1034,9 +1057,6 @@ int colvarmodule::analyze() 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 = variables_active()->begin(); cvi != variables_active()->end(); @@ -1089,7 +1109,6 @@ int colvarmodule::end_of_step() int colvarmodule::setup() { if (this->size() == 0) return cvm::get_error(); - // loop over all components of all colvars to reset masses of all groups for (std::vector::iterator cvi = variables()->begin(); cvi != variables()->end(); cvi++) { (*cvi)->setup(); @@ -1248,13 +1267,13 @@ std::istream & colvarmodule::read_restart(std::istream &is) std::string restart_conf; if (is >> colvarparse::read_block("configuration", restart_conf)) { parse->get_keyval(restart_conf, "step", - it_restart, (size_t) 0, - colvarparse::parse_silent); + it_restart, static_cast(0), + colvarparse::parse_restart); it = it_restart; std::string restart_version; parse->get_keyval(restart_conf, "version", restart_version, std::string(""), - colvarparse::parse_silent); + colvarparse::parse_restart); if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) { cvm::log("This state file was generated with version "+restart_version+"\n"); } @@ -1600,14 +1619,16 @@ std::ostream & colvarmodule::write_traj(std::ostream &os) } -void cvm::log(std::string const &message) +void cvm::log(std::string const &message, int min_log_level) { + if (cvm::log_level() < min_log_level) return; // allow logging when the module is not fully initialized size_t const d = (cvm::main() != NULL) ? depth() : 0; - if (d > 0) + if (d > 0) { proxy->log((std::string(2*d, ' '))+message); - else + } else { proxy->log(message); + } } @@ -1915,14 +1936,190 @@ int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) +template std::string _to_str(T const &x, + size_t width, size_t 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 _to_str_vector(std::vector const &x, + size_t width, size_t 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(); +} + + + +std::string colvarmodule::to_str(std::string const &x) +{ + return std::string("\"")+x+std::string("\""); +} + +std::string colvarmodule::to_str(char const *x) +{ + return std::string("\"")+std::string(x)+std::string("\""); +} + +std::string colvarmodule::to_str(bool x) +{ + return (x ? "on" : "off"); +} + +std::string colvarmodule::to_str(int const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(size_t const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(long int const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(step_number const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(cvm::real const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(cvm::rvector const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(cvm::quaternion const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(colvarvalue const &x, + size_t width, size_t prec) +{ + return _to_str(x, width, prec); +} + +std::string colvarmodule::to_str(cvm::vector1d const &x, + size_t width, size_t prec) +{ + return _to_str< cvm::vector1d >(x, width, prec); +} + +std::string colvarmodule::to_str(cvm::matrix2d const &x, + size_t width, size_t prec) +{ + return _to_str< cvm::matrix2d >(x, width, prec); +} + + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + +std::string colvarmodule::to_str(std::vector const &x, + size_t width, size_t prec) +{ + return _to_str_vector(x, width, prec); +} + + +std::string cvm::wrap_string(std::string const &s, size_t nchars) +{ + if (!s.size()) { + return std::string(nchars, ' '); + } else { + return ( (s.size() <= nchars) ? + (s+std::string(nchars-s.size(), ' ')) : + (std::string(s, 0, nchars)) ); + } +} + + // shared pointer to the proxy object colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07; int colvarmodule::errorCode = 0; -long colvarmodule::it = 0; -long colvarmodule::it_restart = 0; +int colvarmodule::log_level_ = 10; +cvm::step_number colvarmodule::it = 0; +cvm::step_number colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; bool colvarmodule::use_scripted_forces = false; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 99b797627e..30620a2527 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -10,6 +10,8 @@ #ifndef COLVARMODULE_H #define COLVARMODULE_H +#include + #include "colvars_version.h" #ifndef COLVARS_DEBUG @@ -55,6 +57,7 @@ class colvar; class colvarbias; class colvarproxy; class colvarscript; +class colvarvalue; /// \brief Collective variables module (main class) @@ -88,10 +91,16 @@ public: // TODO colvarscript should be unaware of colvarmodule's internals friend class colvarscript; + /// Use a 64-bit integer to store the step number + typedef long long step_number; + /// Defining an abstract real number allows to switch precision typedef double real; - /// Override std::pow with a product for n integer + + // Math functions + + /// Override the STL pow() with a product for n integer static inline real integer_power(real const &x, int const n) { // Original code: math_special.h in LAMMPS @@ -105,8 +114,68 @@ public: return (n > 0) ? yy : 1.0/yy; } - /// Residue identifier - typedef int residue_id; + /// Reimplemented to work around MS compiler issues + static inline real pow(real const &x, real const &y) + { + return ::pow(static_cast(x), static_cast(y)); + } + + /// Reimplemented to work around MS compiler issues + static inline real floor(real const &x) + { + return ::floor(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real fabs(real const &x) + { + return ::fabs(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real sqrt(real const &x) + { + return ::sqrt(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real sin(real const &x) + { + return ::sin(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real cos(real const &x) + { + return ::cos(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real acos(real const &x) + { + return ::acos(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues + static inline real atan2(real const &x, real const &y) + { + return ::atan2(static_cast(x), static_cast(y)); + } + + /// Reimplemented to work around MS compiler issues + static inline real exp(real const &x) + { + return ::exp(static_cast(x)); + } + + /// Reimplemented to work around MS compiler issues. Note: log() is + /// currently defined as the text logging function, but this can be changed + /// at a later time + static inline real logn(real const &x) + { + return ::log(static_cast(x)); + } + class rvector; template class vector1d; @@ -114,6 +183,10 @@ public: class quaternion; class rotation; + + /// Residue identifier + typedef int residue_id; + /// \brief Atom position (different type name from rvector, to make /// possible future PBC-transparent implementations) typedef rvector atom_pos; @@ -150,19 +223,19 @@ public: /// Current step number - static long it; + static step_number it; /// Starting step number for this run - static long it_restart; + static step_number it_restart; /// Return the current step number from the beginning of this run - static inline long step_relative() + static inline step_number step_relative() { return it - it_restart; } /// Return the current step number from the beginning of the whole /// calculation - static inline long step_absolute() + static inline step_number step_absolute() { return it; } @@ -203,9 +276,10 @@ private: std::vector named_atom_groups; public: /// Register a named atom group into named_atom_groups - inline void register_named_atom_group(atom_group * ag) { - named_atom_groups.push_back(ag); - } + void register_named_atom_group(atom_group *ag); + + /// Remove a named atom group from named_atom_groups + void unregister_named_atom_group(atom_group *ag); /// Array of collective variables std::vector *variables(); @@ -254,8 +328,7 @@ public: /// \brief How many objects are configured yet? size_t size() const; - /// \brief Constructor \param config_name Configuration file name - /// \param restart_name (optional) Restart file name + /// \brief Constructor colvarmodule(colvarproxy *proxy); /// Destructor @@ -265,6 +338,7 @@ public: int reset(); /// Open a config file, load its contents, and pass it to config_string() + /// \param config_file_name Configuration file name int read_config_file(char const *config_file_name); /// \brief Parse a config string assuming it is a complete configuration @@ -431,26 +505,92 @@ public: long traj_read_begin, long traj_read_end); - /// 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); + /// Convert to string for output purposes + static std::string to_str(char const *s); + + /// Convert to string for output purposes + static std::string to_str(std::string const &s); + + /// Convert to string for output purposes + static std::string to_str(bool x); + + /// Convert to string for output purposes + static std::string to_str(int const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(size_t const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(long int const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(step_number const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(real const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(rvector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(quaternion const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(colvarvalue const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(vector1d const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(matrix2d const &x, + size_t width = 0, size_t prec = 0); + + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t prec = 0); + + /// Convert to string for output purposes + static std::string to_str(std::vector const &x, + size_t width = 0, size_t 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)) ); - } + static std::string wrap_string(std::string const &s, + size_t nchars); /// Number of characters to represent a time step static size_t const it_width; @@ -485,7 +625,9 @@ public: static void request_total_force(); /// Print a message to the main log - static void log(std::string const &message); + /// \param message Message to print + /// \param min_log_level Only print if cvm::log_level() >= min_log_level + static void log(std::string const &message, int min_log_level = 10); /// Print a message to the main log and exit with error code static int fatal_error(std::string const &message); @@ -493,6 +635,50 @@ public: /// Print a message to the main log and set global error code static int error(std::string const &message, int code = COLVARS_ERROR); +private: + + /// Level of logging requested by the user + static int log_level_; + +public: + + /// Level of logging requested by the user + static inline int log_level() + { + return log_level_; + } + + /// Level at which initialization messages are logged + static inline int log_init_messages() + { + return 1; + } + + /// Level at which a keyword's user-provided value is logged + static inline int log_user_params() + { + return 2; + } + + /// Level at which a keyword's default value is logged + static inline int log_default_params() + { + return 3; + } + + /// Level at which output-file operations are logged + static inline int log_output_files() + { + return 4; + } + + /// Level at which input-file operations (configuration, state) are logged + static inline int log_input_files() + { + return 5; + } + + // Replica exchange commands. static bool replica_enabled(); static int replica_index(); @@ -630,41 +816,4 @@ 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(); -} - - #endif diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index d8b3a359cd..10dce911ee 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -10,6 +10,7 @@ #include #include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -20,15 +21,21 @@ char const * const colvarparse::white_space = " \t"; -// definition of single-value keyword parsers +namespace { -template bool colvarparse::_get_keyval_scalar_(std::string const &conf, - char const *key, - TYPE &value, - TYPE const &def_value, - Parse_Mode const parse_mode) + // Avoid having to put the bool assignment in the template :-( + void set_bool(void *p, bool val) + { + bool *v = reinterpret_cast(p); + *v = val; + } + +} + + +bool colvarparse::get_key_string_value(std::string const &conf, + char const *key, std::string &data) { - std::string data; bool b_found = false, b_found_any = false; size_t save_pos = 0, found_count = 0; @@ -48,29 +55,163 @@ template bool colvarparse::_get_keyval_scalar_(std::string const std::string(key)+"\".\n", INPUT_ERROR); } + return b_found_any; +} + + +template +void colvarparse::mark_key_set_user(std::string const &key_str, + TYPE const &value, + Parse_Mode const &parse_mode) +{ + key_set_modes[to_lower_cppstr(key_str)] = key_set_user; + if (parse_mode & parse_echo) { + cvm::log("# "+key_str+" = "+cvm::to_str(value)+"\n", + cvm::log_user_params()); + } +} + + +template +void colvarparse::mark_key_set_default(std::string const &key_str, + TYPE const &def_value, + Parse_Mode const &parse_mode) +{ + key_set_modes[to_lower_cppstr(key_str)] = key_set_default; + if (parse_mode & parse_echo_default) { + cvm::log("# "+key_str+" = "+cvm::to_str(def_value)+ + " [default]\n", cvm::log_default_params()); + } +} + + +void colvarparse::error_key_required(std::string const &key_str, + Parse_Mode const &parse_mode) +{ + if (key_already_set(key_str)) { + return; + } + if (parse_mode & parse_restart) { + cvm::error("Error: keyword \""+key_str+ + "\" is missing from the restart.\n", INPUT_ERROR); + } else { + cvm::error("Error: keyword \""+key_str+ + "\" is required.\n", INPUT_ERROR); + } +} + + +template +int colvarparse::_get_keyval_scalar_value_(std::string const &key_str, + std::string const &data, + TYPE &value, + TYPE const &def_value) +{ + std::istringstream is(data); + size_t value_count = 0; + TYPE x(def_value); + + while (is >> x) { + value = x; + value_count++; + } + + if (value_count == 0) { + return cvm::error("Error: in parsing \""+ + key_str+"\".\n", INPUT_ERROR); + } + + if (value_count > 1) { + return cvm::error("Error: multiple values " + "are not allowed for keyword \""+ + key_str+"\".\n", INPUT_ERROR); + } + + return COLVARS_OK; +} + + +template<> +int colvarparse::_get_keyval_scalar_value_(std::string const &key_str, + std::string const &data, + bool &value, + bool const &def_value) +{ + if ( (data == std::string("on")) || + (data == std::string("yes")) || + (data == std::string("true")) ) { + set_bool(reinterpret_cast(&value), true); + } else if ( (data == std::string("off")) || + (data == std::string("no")) || + (data == std::string("false")) ) { + set_bool(reinterpret_cast(&value), false); + } else { + return cvm::error("Error: boolean values only are allowed " + "for \""+key_str+"\".\n", INPUT_ERROR); + } + return COLVARS_OK; +} + + +template +int colvarparse::_get_keyval_scalar_novalue_(std::string const &key_str, + TYPE &value, + Parse_Mode const &parse_mode) +{ + return cvm::error("Error: improper or missing value " + "for \""+key_str+"\".\n", INPUT_ERROR); +} + +template<> +int colvarparse::_get_keyval_scalar_novalue_(std::string const &key_str, + bool &value, + Parse_Mode const &parse_mode) +{ + set_bool(reinterpret_cast(&value), true); + mark_key_set_user(key_str, value, parse_mode); + return COLVARS_OK; +} + + +template +bool colvarparse::_get_keyval_scalar_(std::string const &conf, + char const *key, + TYPE &value, + TYPE const &def_value, + Parse_Mode const &parse_mode) +{ + std::string const key_str(key); + + std::string data; + bool const b_found_any = get_key_string_value(conf, key, data); + if (data.size()) { - std::istringstream is(data); - TYPE x(def_value); - if (is >> x) { - value = x; + + _get_keyval_scalar_value_(key_str, data, value, def_value); + + mark_key_set_user(key_str, value, parse_mode); + + } else { // No string value + + if (b_found_any) { + + _get_keyval_scalar_novalue_(key_str, value, parse_mode); + } else { - cvm::error("Error: in parsing \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - cvm::to_str(value)+"\n"); - } - } else { - if (b_found_any) { - cvm::error("Error: improper or missing value " - "for \""+std::string(key)+"\".\n", INPUT_ERROR); - } - value = def_value; - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - cvm::to_str(def_value)+" [default]\n"); + if (parse_mode & parse_required) { + if (cvm::debug()) { + cvm::log("get_keyval, parse_required = "+cvm::to_str(parse_mode & parse_required)+ + "\n"); + } + error_key_required(key_str, parse_mode); + return false; + } + + if ( (parse_mode & parse_override) || !(key_already_set(key)) ) { + value = def_value; + mark_key_set_default(key_str, value, parse_mode); + } } } @@ -78,96 +219,17 @@ template bool colvarparse::_get_keyval_scalar_(std::string const } -bool colvarparse::_get_keyval_scalar_string_(std::string const &conf, - char const *key, - std::string &value, - std::string const &def_value, - Parse_Mode const parse_mode) +template +bool colvarparse::_get_keyval_vector_(std::string const &conf, + char const *key, + std::vector &values, + std::vector const &def_values, + Parse_Mode const &parse_mode) { + std::string const key_str(key); + std::string data; - bool b_found = false, b_found_any = false; - size_t save_pos = 0, found_count = 0; - - do { - std::string data_this = ""; - b_found = key_lookup(conf, key, &data_this, &save_pos); - if (b_found) { - if (!b_found_any) - b_found_any = true; - found_count++; - data = data_this; - } - } while (b_found); - - if (found_count > 1) { - cvm::error("Error: found more than one instance of \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } - - if (data.size()) { - std::istringstream is(data); - size_t data_count = 0; - std::string x(def_value); - while (is >> x) { - value = x; - data_count++; - } - if (data_count == 0) - cvm::error("Error: in parsing \""+ - std::string(key)+"\".\n", INPUT_ERROR); - if (data_count > 1) { - cvm::error("Error: multiple values " - "are not allowed for keyword \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = \""+ - cvm::to_str(value)+"\"\n"); - } - } else { - - if (b_found_any) { - cvm::error("Error: improper or missing value " - "for \""+std::string(key)+"\".\n", INPUT_ERROR); - } - value = def_value; - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = \""+ - cvm::to_str(def_value)+"\" [default]\n"); - } - } - - return b_found_any; -} - - -// multiple-value keyword parsers - -template bool colvarparse::_get_keyval_vector_(std::string const &conf, - char const *key, - std::vector &values, - std::vector const &def_values, - Parse_Mode const parse_mode) -{ - std::string data; - bool b_found = false, b_found_any = false; - size_t save_pos = 0, found_count = 0; - - do { - std::string data_this = ""; - b_found = key_lookup(conf, key, &data_this, &save_pos); - if (b_found) { - if (!b_found_any) - b_found_any = true; - found_count++; - data = data_this; - } - } while (b_found); - - if (found_count > 1) { - cvm::error("Error: found more than one instance of \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } + bool const b_found_any = get_key_string_value(conf, key, data); if (data.size()) { std::istringstream is(data); @@ -175,10 +237,11 @@ template bool colvarparse::_get_keyval_vector_(std::string const if (values.size() == 0) { std::vector x; - if (def_values.size()) + if (def_values.size()) { x = def_values; - else + } else { x.assign(1, TYPE()); + } for (size_t i = 0; ( is >> x[ ((i bool colvarparse::_get_keyval_vector_(std::string const values[i] = x; } else { cvm::error("Error: in parsing \""+ - std::string(key)+"\".\n", INPUT_ERROR); + key_str+"\".\n", INPUT_ERROR); } } } - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - cvm::to_str(values)+"\n"); - } + mark_key_set_user< std::vector >(key_str, values, parse_mode); } else { if (b_found_any) { cvm::error("Error: improper or missing values for \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } + key_str+"\".\n", INPUT_ERROR); + } else { - for (size_t i = 0; i < values.size(); i++) - values[i] = def_values[ (i > def_values.size()) ? 0 : i ]; + if ((values.size() > 0) && (values.size() != def_values.size())) { + cvm::error("Error: the number of default values for \""+ + key_str+"\" is different from the number of " + "current values.\n", BUG_ERROR); + } + + if (parse_mode & parse_required) { + error_key_required(key_str, parse_mode); + return false; + } + + if ( (parse_mode & parse_override) || !(key_already_set(key)) ) { + for (size_t i = 0; i < values.size(); i++) { + values[i] = def_values[i]; + } + mark_key_set_default< std::vector >(key_str, def_values, + parse_mode); + } - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - cvm::to_str(def_values)+" [default]\n"); } } @@ -255,13 +328,22 @@ bool colvarparse::get_keyval(std::string const &conf, return _get_keyval_scalar_(conf, key, value, def_value, parse_mode); } +bool colvarparse::get_keyval(std::string const &conf, + char const *key, + cvm::step_number &value, + cvm::step_number const &def_value, + Parse_Mode const parse_mode) +{ + return _get_keyval_scalar_(conf, key, value, def_value, parse_mode); +} + bool colvarparse::get_keyval(std::string const &conf, char const *key, std::string &value, std::string const &def_value, Parse_Mode const parse_mode) { - return _get_keyval_scalar_string_(conf, key, value, def_value, parse_mode); + return _get_keyval_scalar_(conf, key, value, def_value, parse_mode); } bool colvarparse::get_keyval(std::string const &conf, @@ -300,66 +382,13 @@ bool colvarparse::get_keyval(std::string const &conf, return _get_keyval_scalar_(conf, key, value, def_value, parse_mode); } - bool colvarparse::get_keyval(std::string const &conf, char const *key, bool &value, bool const &def_value, Parse_Mode const parse_mode) { - std::string data; - bool b_found = false, b_found_any = false; - size_t save_pos = 0, found_count = 0; - - do { - std::string data_this = ""; - b_found = key_lookup(conf, key, &data_this, &save_pos); - if (b_found) { - if (!b_found_any) - b_found_any = true; - found_count++; - data = data_this; - } - } while (b_found); - - if (found_count > 1) { - cvm::error("Error: found more than one instance of \""+ - std::string(key)+"\".\n", INPUT_ERROR); - } - - if (data.size()) { - if ( (data == std::string("on")) || - (data == std::string("yes")) || - (data == std::string("true")) ) { - value = true; - } else if ( (data == std::string("off")) || - (data == std::string("no")) || - (data == std::string("false")) ) { - value = false; - } else - cvm::error("Error: boolean values only are allowed " - "for \""+std::string(key)+"\".\n", INPUT_ERROR); - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - (value ? "on" : "off")+"\n"); - } - } else { - - if (b_found_any) { - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = on\n"); - } - value = true; - } else { - value = def_value; - if (parse_mode != parse_silent) { - cvm::log("# "+std::string(key)+" = "+ - (def_value ? "on" : "off")+" [default]\n"); - } - } - } - - return b_found_any; + return _get_keyval_scalar_(conf, key, value, def_value, parse_mode); } @@ -440,15 +469,27 @@ bool colvarparse::get_keyval(std::string const &conf, void colvarparse::add_keyword(char const *key) { - for (std::list::iterator ki = allowed_keywords.begin(); - ki != allowed_keywords.end(); ki++) { - if (to_lower_cppstr(std::string(key)) == *ki) - return; + std::string const key_str_lower(to_lower_cppstr(std::string(key))); + + if (key_set_modes.find(key_str_lower) != key_set_modes.end()) { + return; } - // not found in the list - // if (cvm::debug()) - // cvm::log("Registering a new keyword, \""+std::string (key)+"\".\n"); - allowed_keywords.push_back(to_lower_cppstr(std::string(key))); + + key_set_modes[key_str_lower] = key_not_set; + + allowed_keywords.push_back(key_str_lower); +} + + +bool colvarparse::key_already_set(std::string const &key_str) +{ + std::string const key_str_lower(to_lower_cppstr(key_str)); + + if (key_set_modes.find(key_str_lower) == key_set_modes.end()) { + return false; + } + + return (key_set_modes[key_str_lower] > 0); } @@ -457,6 +498,10 @@ void colvarparse::strip_values(std::string &conf) size_t offset = 0; data_begin_pos.sort(); data_end_pos.sort(); + std::list::iterator data_begin_pos_last = std::unique(data_begin_pos.begin(), data_begin_pos.end()); + data_begin_pos.erase(data_begin_pos_last, data_begin_pos.end()); + std::list::iterator data_end_pos_last = std::unique(data_end_pos.begin(), data_end_pos.end()); + data_end_pos.erase(data_end_pos_last, data_end_pos.end()); std::list::iterator data_begin = data_begin_pos.begin(); std::list::iterator data_end = data_end_pos.begin(); @@ -464,28 +509,16 @@ void colvarparse::strip_values(std::string &conf) for ( ; (data_begin != data_begin_pos.end()) && (data_end != data_end_pos.end()) ; data_begin++, data_end++) { - - // std::cerr << "data_begin, data_end " - // << *data_begin << ", " << *data_end - // << "\n"; - size_t const nchars = *data_end-*data_begin; - - // std::cerr << "conf[data_begin:data_end] = \"" - // << std::string (conf, *data_begin - offset, nchars) - // << "\"\n"; - conf.erase(*data_begin - offset, nchars); offset += nchars; - - // std::cerr << ("Stripped config = \"\n"+conf+"\"\n"); - } } void colvarparse::clear_keyword_registry() { + key_set_modes.clear(); allowed_keywords.clear(); data_begin_pos.clear(); data_end_pos.clear(); diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index 28ad3c052b..8501ee8c14 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -12,6 +12,7 @@ #include #include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -24,35 +25,9 @@ /// need to parse input inherit from this class colvarparse { -protected: - - /// \brief List of legal keywords for this object: this is updated - /// by each call to colvarparse::get_keyval() or - /// colvarparse::key_lookup() - std::list allowed_keywords; - - /// \brief List of delimiters for the values of each keyword in the - /// configuration string; all keywords will be stripped of their - /// values before the keyword check is performed - std::list data_begin_pos; - - /// \brief List of delimiters for the values of each keyword in the - /// configuration string; all keywords will be stripped of their - /// values before the keyword check is performed - std::list data_end_pos; - - /// \brief Add a new valid keyword to the list - void add_keyword(char const *key); - - /// \brief Remove all the values from the config string - void strip_values(std::string &conf); - - /// \brief Configuration string of the object (includes comments) - std::string config_string; - public: - + /// Default constructor inline colvarparse() { init(); @@ -88,14 +63,23 @@ public: /// How a keyword is parsed in a string enum Parse_Mode { - /// \brief(default) Read the first instance of a keyword (if - /// any), report its value, and print a warning when there is more - /// than one - parse_normal, - /// \brief Like parse_normal, but don't send any message to the log - /// (useful e.g. in restart files when such messages are very - /// numerous and redundant) - parse_silent + /// Zero for all flags + parse_null = 0, + /// Print the value of a keyword if it is given + parse_echo = (1<<1), + /// Print the default value of a keyword, if it is NOT given + parse_echo_default = (1<<2), + /// Do not print the keyword + parse_silent = 0, + /// Raise error if the keyword is not provided + parse_required = (1<<16), + /// Successive calls to get_keyval() will override the previous values + /// when the keyword is not given any more + parse_override = (1<<17), + /// The call is being executed from a read_restart() function + parse_restart = (1<<18), + /// Alias for old default behavior (should be phased out) + parse_normal = (1<<2) | (1<<1) | (1<<17) }; /// \brief Check that all the keywords within "conf" are in the list @@ -146,6 +130,11 @@ public: long &value, long const &def_value = 0, Parse_Mode const parse_mode = parse_normal); + bool get_keyval(std::string const &conf, + char const *key, + cvm::step_number &value, + cvm::step_number const &def_value = 0, + Parse_Mode const parse_mode = parse_normal); bool get_keyval(std::string const &conf, char const *key, std::string &value, @@ -219,23 +208,57 @@ public: protected: - // Templates - template bool _get_keyval_scalar_(std::string const &conf, - char const *key, - TYPE &value, - TYPE const &def_value, - Parse_Mode const parse_mode); - bool _get_keyval_scalar_string_(std::string const &conf, - char const *key, - std::string &value, - std::string const &def_value, - Parse_Mode const parse_mode); + /// Get the string value of a keyword, and save it for later parsing + bool get_key_string_value(std::string const &conf, + char const *key, std::string &data); - template bool _get_keyval_vector_(std::string const &conf, - char const *key, - std::vector &values, - std::vector const &def_values, - Parse_Mode const parse_mode); + /// Template for single-value keyword parsers + template + bool _get_keyval_scalar_(std::string const &conf, + char const *key, + TYPE &value, + TYPE const &def_value, + Parse_Mode const &parse_mode); + + /// Template for multiple-value keyword parsers + template + bool _get_keyval_vector_(std::string const &conf, + char const *key, + std::vector &values, + std::vector const &def_values, + Parse_Mode const &parse_mode); + + /// Extract the value of a variable from a string + template + int _get_keyval_scalar_value_(std::string const &key_str, + std::string const &data, + TYPE &value, + TYPE const &def_value); + + /// Handle the case where the user provides a keyword without value + template + int _get_keyval_scalar_novalue_(std::string const &key_str, + TYPE &value, + Parse_Mode const &parse_mode); + + /// Record that the keyword has just been user-defined + template + void mark_key_set_user(std::string const &key_str, + TYPE const &value, + Parse_Mode const &parse_mode); + + /// Record that the keyword has just been set to its default value + template + void mark_key_set_default(std::string const &key_str, + TYPE const &def_value, + Parse_Mode const &parse_mode); + + /// Raise error condition due to the keyword being required! + void error_key_required(std::string const &key_str, + Parse_Mode const &parse_mode); + + /// True if the keyword has been set already + bool key_already_set(std::string const &key_str); public: @@ -286,7 +309,7 @@ public: size_t *save_pos = NULL); /// \brief Reads a configuration line, adds it to config_string, and returns - /// the stream \param is Input stream \param s String that will hold the + /// the stream \param is Input stream \param line String that will hold the /// configuration line, with comments stripped std::istream & read_config_line(std::istream &is, std::string &line); @@ -299,7 +322,51 @@ public: /// from this position static int check_braces(std::string const &conf, size_t const start_pos); +protected: + + /// \brief List of legal keywords for this object: this is updated + /// by each call to colvarparse::get_keyval() or + /// colvarparse::key_lookup() + std::list allowed_keywords; + + /// How a keyword has been set + enum key_set_mode { + key_not_set = 0, + key_set_user = 1, + key_set_default = 2 + }; + + /// Track which keywords have been already set, and how + std::map key_set_modes; + + /// \brief List of delimiters for the values of each keyword in the + /// configuration string; all keywords will be stripped of their + /// values before the keyword check is performed + std::list data_begin_pos; + + /// \brief List of delimiters for the values of each keyword in the + /// configuration string; all keywords will be stripped of their + /// values before the keyword check is performed + std::list data_end_pos; + + /// \brief Add a new valid keyword to the list + void add_keyword(char const *key); + + /// \brief Remove all the values from the config string + void strip_values(std::string &conf); + + /// \brief Configuration string of the object (includes comments) + std::string config_string; + }; +/// Bitwise OR between two Parse_mode flags +inline colvarparse::Parse_Mode operator | (colvarparse::Parse_Mode const &mode1, + colvarparse::Parse_Mode const &mode2) +{ + return static_cast(static_cast(mode1) | + static_cast(mode2)); +} + #endif diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp index da9257eee2..5f8e82d30a 100644 --- a/lib/colvars/colvarproxy.cpp +++ b/lib/colvars/colvarproxy.cpp @@ -60,7 +60,7 @@ bool colvarproxy_system::total_forces_same_step() const inline int round_to_integer(cvm::real x) { - return std::floor(x+0.5); + return cvm::floor(x+0.5); } @@ -129,7 +129,10 @@ cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1, -colvarproxy_atoms::colvarproxy_atoms() {} +colvarproxy_atoms::colvarproxy_atoms() +{ + updated_masses_ = updated_charges_ = false; +} colvarproxy_atoms::~colvarproxy_atoms() @@ -544,7 +547,7 @@ int colvarproxy_script::run_colvar_gradient_callback( colvarproxy_tcl::colvarproxy_tcl() { - _tcl_interp = NULL; + tcl_interp_ = NULL; } @@ -573,7 +576,7 @@ char const *colvarproxy_tcl::tcl_obj_to_str(unsigned char *obj) int colvarproxy_tcl::tcl_run_force_callback() { #if defined(COLVARS_TCL) - Tcl_Interp *const tcl_interp = reinterpret_cast(_tcl_interp); + Tcl_Interp *const tcl_interp = reinterpret_cast(tcl_interp_); std::string cmd = std::string("calc_colvar_forces ") + cvm::to_str(cvm::step_absolute()); int err = Tcl_Eval(tcl_interp, cmd.c_str()); @@ -596,7 +599,7 @@ int colvarproxy_tcl::tcl_run_colvar_callback( { #if defined(COLVARS_TCL) - Tcl_Interp *const tcl_interp = reinterpret_cast(_tcl_interp); + Tcl_Interp *const tcl_interp = reinterpret_cast(tcl_interp_); size_t i; std::string cmd = std::string("calc_") + name; for (i = 0; i < cvc_values.size(); i++) { @@ -633,7 +636,7 @@ int colvarproxy_tcl::tcl_run_colvar_gradient_callback( { #if defined(COLVARS_TCL) - Tcl_Interp *const tcl_interp = reinterpret_cast(_tcl_interp); + Tcl_Interp *const tcl_interp = reinterpret_cast(tcl_interp_); size_t i; std::string cmd = std::string("calc_") + name + "_gradient"; for (i = 0; i < cvc_values.size(); i++) { diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 3ff5f3c697..3bbdfe522e 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -29,7 +29,7 @@ /// /// To interface to a new MD engine, the simplest solution is to derive a new /// class from \link colvarproxy \endlink. Currently implemented are: \link -/// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link +/// colvarproxy_lammps \endlink, \link colvarproxy_namd \endlink, \link /// colvarproxy_vmd \endlink. @@ -227,11 +227,15 @@ public: inline std::vector *modify_atom_masses() { + // assume that we are requesting masses to change them + updated_masses_ = true; return &atoms_masses; } inline std::vector *modify_atom_charges() { + // assume that we are requesting charges to change them + updated_charges_ = true; return &atoms_charges; } @@ -250,6 +254,18 @@ public: return &atoms_new_colvar_forces; } + /// Record whether masses have been updated + inline bool updated_masses() const + { + return updated_masses_; + } + + /// Record whether masses have been updated + inline bool updated_charges() const + { + return updated_charges_; + } + protected: /// \brief Array of 0-based integers used to uniquely associate atoms @@ -268,6 +284,9 @@ protected: /// \brief Forces applied from colvars, to be communicated to the MD integrator std::vector atoms_new_colvar_forces; + /// Whether the masses and charges have been updated from the host code + bool updated_masses_, updated_charges_; + /// Used by all init_atom() functions: create a slot for an atom not /// requested yet; returns the index in the arrays int add_atom_slot(int atom_id); @@ -522,7 +541,7 @@ public: protected: /// Pointer to Tcl interpreter object - void *_tcl_interp; + void *tcl_interp_; /// Set Tcl pointers virtual void init_tcl_pointers(); diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index bd84d077d7..2521fdf872 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,5 +1,5 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2018-11-16" +#define COLVARS_VERSION "2019-04-26" // This file is part of the Collective Variables module (Colvars). // The original version of Colvars and its updates are located at: // https://github.com/colvars/colvars diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index c9fe0497a9..1ad3283337 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -137,6 +137,10 @@ int colvarscript::run(int objc, unsigned char *const objv[]) if (cmd == "update") { error_code |= proxy->update_input(); + if (error_code) { + result += "Error updating the Colvars module.\n"; + return error_code; + } error_code |= colvars->calc(); error_code |= proxy->update_output(); if (error_code) { @@ -273,6 +277,10 @@ int colvarscript::run(int objc, unsigned char *const objv[]) int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[]) { + if (objc < 3) { + result = "Missing arguments"; + return COLVARSCRIPT_ERROR; + } std::string const subcmd(obj_to_str(objv[2])); if (subcmd == "value") { @@ -323,6 +331,47 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[]) return COLVARS_OK; } + if (subcmd == "getatomgroups") { + std::vector > lists = cv->get_atom_lists(); + std::vector >::iterator li = lists.begin(); + + for ( ; li != lists.end(); ++li) { + result += "{"; + std::vector::iterator lj = (*li).begin(); + for ( ; lj != (*li).end(); ++lj) { + result += cvm::to_str(*lj); + result += " "; + } + result += "} "; + } + return COLVARS_OK; + } + + if (subcmd == "getatomids") { + std::vector::iterator li = cv->atom_ids.begin(); + + for ( ; li != cv->atom_ids.end(); ++li) { + result += cvm::to_str(*li); + result += " "; + } + return COLVARS_OK; + } + + if (subcmd == "getgradients") { + std::vector::iterator li = cv->atomic_gradients.begin(); + + for ( ; li != cv->atomic_gradients.end(); ++li) { + result += "{"; + int j; + for (j = 0; j < 3; ++j) { + result += cvm::to_str((*li)[j]); + result += " "; + } + result += "} "; + } + return COLVARS_OK; + } + if (subcmd == "getappliedforce") { result = (cv->applied_force()).to_simple_string(); return COLVARS_OK; @@ -410,6 +459,10 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[]) int colvarscript::proc_bias(colvarbias *b, int objc, unsigned char *const objv[]) { + if (objc < 3) { + result = "Missing arguments"; + return COLVARSCRIPT_ERROR; + } std::string const subcmd(obj_to_str(objv[2])); if (subcmd == "energy") { diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index 313dbd661b..341cb1f72c 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -191,8 +191,9 @@ inline static colvarbias *colvarbias_obj(void *pobj) #ifdef COLVARSCRIPT_CPP #define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \ - int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \ - int objc, unsigned char *const objv[]) \ + extern "C" int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \ + int objc, \ + unsigned char *const objv[]) \ { \ colvarscript *script = colvarscript_obj(); \ script->clear_results(); \ diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 2b45d77e07..b81a943eab 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -19,6 +19,8 @@ bool colvarmodule::rotation::monitor_crossings = false; cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-02; +namespace { + /// Numerical recipes diagonalization static int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot); @@ -28,6 +30,7 @@ static int eigsrt(cvm::real *d, cvm::real **v); /// Transpose the matrix static int transpose(cvm::real **v); +} std::string cvm::rvector::to_simple_string() const @@ -245,13 +248,11 @@ cvm::quaternion::position_derivative_inner(cvm::rvector const &pos, // 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, - cvm::matrix2d &S) +void colvarmodule::rotation::build_correlation_matrix( + std::vector const &pos1, + std::vector const &pos2) { // build the correlation matrix - C.resize(3, 3); - C.reset(); size_t i; for (i = 0; i < pos1.size(); i++) { C.xx() += pos1[i].x * pos2[i].x; @@ -264,7 +265,11 @@ void colvarmodule::rotation::build_matrix(std::vector const &pos1 C.zy() += pos1[i].z * pos2[i].y; C.zz() += pos1[i].z * pos2[i].z; } +} + +void colvarmodule::rotation::compute_overlap_matrix() +{ // 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(); @@ -286,37 +291,38 @@ void colvarmodule::rotation::build_matrix(std::vector const &pos1 } -void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d &S, - cvm::vector1d &S_eigval, - cvm::matrix2d &S_eigvec) +void colvarmodule::rotation::diagonalize_matrix( + cvm::matrix2d &m, + cvm::vector1d &eigval, + cvm::matrix2d &eigvec) { - S_eigval.resize(4); - S_eigval.reset(); - S_eigvec.resize(4,4); - S_eigvec.reset(); + eigval.resize(4); + eigval.reset(); + eigvec.resize(4, 4); + eigvec.reset(); // diagonalize int jac_nrot = 0; - if (jacobi(S.c_array(), S_eigval.c_array(), S_eigvec.c_array(), &jac_nrot) != + if (jacobi(m.c_array(), eigval.c_array(), eigvec.c_array(), &jac_nrot) != COLVARS_OK) { cvm::error("Too many iterations in routine jacobi.\n" "This is usually the result of an ill-defined set of atoms for " "rotational alignment (RMSD, rotateReference, etc).\n"); } - eigsrt(S_eigval.c_array(), S_eigvec.c_array()); + eigsrt(eigval.c_array(), eigvec.c_array()); // jacobi saves eigenvectors by columns - transpose(S_eigvec.c_array()); + transpose(eigvec.c_array()); // normalize eigenvectors for (size_t ie = 0; ie < 4; ie++) { cvm::real norm2 = 0.0; size_t i; for (i = 0; i < 4; i++) { - norm2 += S_eigvec[ie][i] * S_eigvec[ie][i]; + norm2 += eigvec[ie][i] * eigvec[ie][i]; } - cvm::real const norm = std::sqrt(norm2); + cvm::real const norm = cvm::sqrt(norm2); for (i = 0; i < 4; i++) { - S_eigvec[ie][i] /= norm; + eigvec[ie][i] /= norm; } } } @@ -324,23 +330,26 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d &S, // Calculate the rotation, plus its derivatives -void colvarmodule::rotation::calc_optimal_rotation(std::vector const &pos1, - std::vector const &pos2) +void colvarmodule::rotation::calc_optimal_rotation( + std::vector const &pos1, + std::vector const &pos2) { - S.resize(4,4); + C.resize(3, 3); + C.reset(); + build_correlation_matrix(pos1, pos2); + + S.resize(4, 4); S.reset(); + compute_overlap_matrix(); - build_matrix(pos1, pos2, S); - - S_backup.resize(4,4); + S_backup.resize(4, 4); S_backup = S; if (b_debug_gradients) { - cvm::log("S = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n"); + cvm::log("S = "+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]; @@ -524,7 +533,7 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector co 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)+ + cvm::to_str(cvm::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"); @@ -544,6 +553,9 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector co #define n 4 + +namespace { + int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot) { int j,iq,ip,i; @@ -567,7 +579,7 @@ int jacobi(cvm::real **a, cvm::real *d, cvm::real **v, int *nrot) sm=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])) + g=100.0*cvm::fabs(a[ip][iq]); + if (i > 4 && (cvm::real)(cvm::fabs(d[ip])+g) == (cvm::real)cvm::fabs(d[ip]) + && (cvm::real)(cvm::fabs(d[iq])+g) == (cvm::real)cvm::fabs(d[iq])) a[ip][iq]=0.0; - else if (std::fabs(a[ip][iq]) > tresh) { + else if (cvm::fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; - if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h)) + if ((cvm::real)(cvm::fabs(h)+g) == (cvm::real)cvm::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)); + t=1.0/(cvm::fabs(theta)+cvm::sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } - c=1.0/std::sqrt(1+t*t); + c=1.0/cvm::sqrt(1+t*t); s=t*c; tau=s/(1.0+c); h=t*a[ip][iq]; @@ -663,5 +675,7 @@ int transpose(cvm::real **v) return COLVARS_OK; } +} + #undef n #undef ROTATE diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index 4ef9557dcb..9973d92e98 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -10,7 +10,6 @@ #ifndef COLVARTYPES_H #define COLVARTYPES_H -#include #include #include "colvarmodule.h" @@ -220,7 +219,7 @@ public: inline cvm::real norm() const { - return std::sqrt(this->norm2()); + return cvm::sqrt(this->norm2()); } inline cvm::real sum() const @@ -801,7 +800,7 @@ public: inline cvm::real norm() const { - return std::sqrt(this->norm2()); + return cvm::sqrt(this->norm2()); } inline cvm::rvector unit() const @@ -1008,17 +1007,17 @@ public: cvm::real theta_in, cvm::real psi_in) { - q0 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) + - (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) ); + q0 = ( (cvm::cos(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) + + (cvm::sin(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) ); - q1 = ( (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) - - (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) ); + q1 = ( (cvm::sin(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) - + (cvm::cos(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) ); - q2 = ( (std::cos(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) + - (std::sin(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) ); + q2 = ( (cvm::cos(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) + + (cvm::sin(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) ); - q3 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::sin(psi_in/2.0)) - - (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::cos(psi_in/2.0)) ); + q3 = ( (cvm::cos(phi_in/2.0)) * (cvm::cos(theta_in/2.0)) * (cvm::sin(psi_in/2.0)) - + (cvm::sin(phi_in/2.0)) * (cvm::sin(theta_in/2.0)) * (cvm::cos(psi_in/2.0)) ); } /// \brief Default constructor @@ -1115,7 +1114,7 @@ public: /// Norm of the quaternion inline cvm::real norm() const { - return std::sqrt(this->norm2()); + return cvm::sqrt(this->norm2()); } /// Return the conjugate quaternion @@ -1177,7 +1176,7 @@ public: } /// \brief Provides the quaternion product. \b NOTE: for the inner - /// product use: \code h.inner (q); \endcode + /// product use: `h.inner (q);` friend inline cvm::quaternion operator * (cvm::quaternion const &h, cvm::quaternion const &q) { @@ -1263,7 +1262,7 @@ public: cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 + this->q2*Q2.q2 + this->q3*Q2.q3; - cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 : + cvm::real const omega = cvm::acos( (cos_omega > 1.0) ? 1.0 : ( (cos_omega < -1.0) ? -1.0 : cos_omega) ); // get the minimum distance: x and -x are the same quaternion @@ -1278,11 +1277,11 @@ public: inline cvm::quaternion dist2_grad(cvm::quaternion const &Q2) const { cvm::real const cos_omega = this->q0*Q2.q0 + this->q1*Q2.q1 + this->q2*Q2.q2 + this->q3*Q2.q3; - cvm::real const omega = std::acos( (cos_omega > 1.0) ? 1.0 : + cvm::real const omega = cvm::acos( (cos_omega > 1.0) ? 1.0 : ( (cos_omega < -1.0) ? -1.0 : cos_omega) ); - cvm::real const sin_omega = std::sin(omega); + cvm::real const sin_omega = cvm::sin(omega); - if (std::fabs(sin_omega) < 1.0E-14) { + if (cvm::fabs(sin_omega) < 1.0E-14) { // return a null 4d vector return cvm::quaternion(0.0, 0.0, 0.0, 0.0); } @@ -1338,14 +1337,16 @@ public: /// \brief Perform gradient tests bool b_debug_gradients; - /// \brief Positions to superimpose: the rotation should brings pos1 - /// into pos2 - std::vector pos1, pos2; - + /// Correlation matrix C (3, 3) cvm::rmatrix C; + /// Overlap matrix S (4, 4) cvm::matrix2d S; + + /// Eigenvalues of S cvm::vector1d S_eigval; + + /// Eigenvectors of S cvm::matrix2d S_eigvec; /// Used for debugging gradients @@ -1404,8 +1405,8 @@ public: : b_debug_gradients(false) { cvm::rvector const axis_n = axis.unit(); - cvm::real const sina = std::sin(angle/2.0); - q = cvm::quaternion(std::cos(angle/2.0), + cvm::real const sina = cvm::sin(angle/2.0); + q = cvm::quaternion(cvm::cos(angle/2.0), sina * axis_n.x, sina * axis_n.y, sina * axis_n.z); } @@ -1437,7 +1438,7 @@ public: inline cvm::real spin_angle(cvm::rvector const &axis) const { cvm::rvector const q_vec = q.get_vector(); - cvm::real alpha = (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0); + cvm::real alpha = (180.0/PI) * 2.0 * cvm::atan2(axis * q_vec, q.q0); while (alpha > 180.0) alpha -= 360; while (alpha < -180.0) alpha += 360; return alpha; @@ -1473,9 +1474,9 @@ public: { cvm::rvector const q_vec = q.get_vector(); cvm::real const alpha = - (180.0/PI) * 2.0 * std::atan2(axis * q_vec, q.q0); + (180.0/PI) * 2.0 * cvm::atan2(axis * q_vec, q.q0); - cvm::real const cos_spin_2 = std::cos(alpha * (PI/180.0) * 0.5); + cvm::real const cos_spin_2 = cvm::cos(alpha * (PI/180.0) * 0.5); cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ? (q.q0 / cos_spin_2) : (0.0) ); @@ -1489,7 +1490,7 @@ public: cvm::rvector const q_vec = q.get_vector(); cvm::real const iprod = axis * q_vec; - cvm::real const cos_spin_2 = std::cos(std::atan2(iprod, q.q0)); + cvm::real const cos_spin_2 = cvm::cos(cvm::atan2(iprod, q.q0)); if (q.q0 != 0.0) { @@ -1529,15 +1530,17 @@ protected: /// eigenvalue crossing) cvm::quaternion q_old; - /// Build the overlap matrix S (used by calc_optimal_rotation()) - void build_matrix(std::vector const &pos1, - std::vector const &pos2, - cvm::matrix2d &S); + /// Build the correlation matrix C (used by calc_optimal_rotation()) + void build_correlation_matrix(std::vector const &pos1, + std::vector const &pos2); - /// Diagonalize the overlap matrix S (used by calc_optimal_rotation()) - void diagonalize_matrix(cvm::matrix2d &S, - cvm::vector1d &S_eigval, - cvm::matrix2d &S_eigvec); + /// Compute the overlap matrix S (used by calc_optimal_rotation()) + void compute_overlap_matrix(); + + /// Diagonalize a given matrix m (used by calc_optimal_rotation()) + static void diagonalize_matrix(cvm::matrix2d &m, + cvm::vector1d &eigval, + cvm::matrix2d &eigvec); }; diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index 86b99ed3d6..accc5defec 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -144,10 +144,10 @@ void colvarvalue::apply_constraints() case colvarvalue::type_quaternionderiv: break; case colvarvalue::type_unit3vector: - rvector_value /= std::sqrt(rvector_value.norm2()); + rvector_value /= cvm::sqrt(rvector_value.norm2()); break; case colvarvalue::type_quaternion: - quaternion_value /= std::sqrt(quaternion_value.norm2()); + quaternion_value /= cvm::sqrt(quaternion_value.norm2()); break; case colvarvalue::type_vector: if (elem_types.size() > 0) { @@ -579,7 +579,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const cvm::rvector const &v1 = this->rvector_value; cvm::rvector const &v2 = x2.rvector_value; cvm::real const cos_t = v1 * v2; - cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t); + cvm::real const sin_t = cvm::sqrt(1.0 - cos_t*cos_t); return colvarvalue( 2.0 * sin_t * cvm::rvector((-1.0) * sin_t * v2.x + cos_t/sin_t * (v1.x - cos_t*v2.x), @@ -630,7 +630,7 @@ colvarvalue const colvarvalue::interpolate(colvarvalue const &x1, break; case colvarvalue::type_unit3vector: case colvarvalue::type_quaternion: - if (interp.norm()/std::sqrt(d2) < 1.0e-6) { + if (interp.norm()/cvm::sqrt(d2) < 1.0e-6) { cvm::error("Error: interpolation between "+cvm::to_str(x1)+" and "+ cvm::to_str(x2)+" with lambda = "+cvm::to_str(lambda)+ " is undefined: result = "+cvm::to_str(interp)+"\n", diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index 25255e2f7c..29b535a1a8 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -17,22 +17,22 @@ /// \brief Value of a collective variable: this is a metatype which /// can be set at runtime. By default it is set to be a scalar /// number, and can be treated as such in all operations (this is -/// done by most \link cvc \endlink implementations). +/// done by most \link colvar::cvc \endlink implementations). /// /// \link colvarvalue \endlink allows \link colvar \endlink to be /// treat different data types. By default, a \link colvarvalue /// \endlink variable is a scalar number. To use it as /// another type, declare and initialize it as -/// \code colvarvalue x(colvarvalue::type_xxx)\endcode, use \link x.type -/// (colvarvalue::type_xxx) \endlink at a later stage, or if unset, -/// assign the type with \code x = y; \endcode, provided y is correctly set. +/// `colvarvalue x(colvarvalue::type_xxx)`, use `x.type (colvarvalue::type_xxx)` +/// at a later stage, or if unset, +/// assign the type with `x = y;`, provided y is correctly set. /// /// All operators (either unary or binary) on a \link /// colvarvalue \endlink object performs one or more checks on the /// \link Type \endlink, except when reading from a stream, when there is no way to -/// detect the \link Type \endlink. To use \code is >> x; \endcode x \b MUST +/// detect the \link Type \endlink. To use `is >> x;` x \b MUST /// already have a type correcly set up for properly parsing the -/// stream. No problem of course with the output streams: \code os << x; \endcode +/// stream. No problem of course with the output streams: `os << x;` /// /// \em Note \em on \em performance: to avoid type checks in a long array of \link /// colvarvalue \endlink objects, use one of the existing "_opt" functions or implement a new one @@ -159,7 +159,7 @@ public: /// \brief If the variable has constraints (e.g. unitvector or /// quaternion), transform it to satisfy them; this function needs /// to be called only when the \link colvarvalue \endlink - /// is calculated outside of \link cvc \endlink objects + /// is calculated outside of \link colvar::cvc \endlink objects void apply_constraints(); /// Get the current type @@ -184,7 +184,7 @@ public: /// Norm of this colvarvalue inline cvm::real norm() const { - return std::sqrt(this->norm2()); + return cvm::sqrt(this->norm2()); } /// Sum of the components of this colvarvalue (if more than one dimension) @@ -728,7 +728,7 @@ inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const case colvarvalue::type_unit3vector: case colvarvalue::type_unit3vectorderiv: // angle between (*this) and x2 is the distance - return std::acos(this->rvector_value * x2.rvector_value) * std::acos(this->rvector_value * x2.rvector_value); + return cvm::acos(this->rvector_value * x2.rvector_value) * cvm::acos(this->rvector_value * x2.rvector_value); case colvarvalue::type_quaternion: case colvarvalue::type_quaternionderiv: // angle between (*this) and x2 is the distance, the quaternion diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index ff74602b41..651999f79c 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -153,18 +153,29 @@ void colvarproxy_lammps::init(const char *conf_file) if (_lmp->update->ntimestep != 0) { cvm::log("Setting initial step number from LAMMPS: "+ cvm::to_str(_lmp->update->ntimestep)+"\n"); - colvars->it = colvars->it_restart = _lmp->update->ntimestep; + colvars->it = colvars->it_restart = + static_cast(_lmp->update->ntimestep); } if (cvm::debug()) { - log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); - log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); - log(cvm::line_marker); - log("Info: done initializing the colvars proxy object.\n"); + cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); + cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); + cvm::log(cvm::line_marker); + cvm::log("Info: done initializing the colvars proxy object.\n"); } } +void colvarproxy_lammps::add_config_file(const char *conf_file) +{ + colvars->read_config_file(conf_file); +} + +void colvarproxy_lammps::add_config_string(const std::string &conf) +{ + colvars->read_config_string(conf); +} + colvarproxy_lammps::~colvarproxy_lammps() { delete _random; @@ -185,7 +196,7 @@ int colvarproxy_lammps::setup() double colvarproxy_lammps::compute() { if (cvm::debug()) { - log(std::string(cvm::line_marker)+ + cvm::log(std::string(cvm::line_marker)+ "colvarproxy_lammps step no. "+ cvm::to_str(_lmp->update->ntimestep)+" [first - last = "+ cvm::to_str(_lmp->update->beginstep)+" - "+ @@ -238,20 +249,20 @@ double colvarproxy_lammps::compute() bias_energy = 0.0; if (cvm::debug()) { - log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); - log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); - log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); + cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); + cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); + cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } // call the collective variable module colvars->calc(); if (cvm::debug()) { - log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); - log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); - log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); + cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); + cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); + cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } return bias_energy; diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index cdd86cbd16..c3d9dbb35f 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -103,6 +103,11 @@ class colvarproxy_lammps : public colvarproxy { // Write files expected from Colvars (called by post_run()) void write_output_files(); + // read additional config from file + void add_config_file(char const *config_filename); + + // read additional config from string + void add_config_string(const std::string &config); // implementation of pure methods from base class public: diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h index edd6778c64..0a4f9fdf4f 100644 --- a/src/USER-COLVARS/colvarproxy_lammps_version.h +++ b/src/USER-COLVARS/colvarproxy_lammps_version.h @@ -1,5 +1,5 @@ #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2018-08-29" +#define COLVARPROXY_VERSION "2019-04-09" // This file is part of the Collective Variables module (Colvars). // The original version of Colvars and its updates are located at: // https://github.com/colvars/colvars diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 545ceb7b0e..ff0e8fb334 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -480,6 +480,31 @@ void FixColvars::one_time_init() /* ---------------------------------------------------------------------- */ +int FixColvars::modify_param(int narg, char **arg) +{ + if (strcmp(arg[0],"configfile") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (me == 0) { + if (! proxy) + error->one(FLERR,"Cannot use fix_modify before initialization"); + proxy->add_config_file(arg[1]); + } + return 2; + } else if (strcmp(arg[0],"config") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (me == 0) { + if (! proxy) + error->one(FLERR,"Cannot use fix_modify before initialization"); + std::string conf(arg[1]); + proxy->add_config_string(conf); + } + return 2; + } + return 0; +} + +/* ---------------------------------------------------------------------- */ + void FixColvars::setup(int vflag) { const tagint * const tag = atom->tag; diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h index 3029ba9db5..a0c197fca4 100644 --- a/src/USER-COLVARS/fix_colvars.h +++ b/src/USER-COLVARS/fix_colvars.h @@ -50,6 +50,7 @@ class FixColvars : public Fix { virtual int setmask(); virtual void init(); virtual void setup(int); + virtual int modify_param(int, char **); virtual void min_setup(int vflag) {setup(vflag);}; virtual void min_post_force(int); virtual void post_force(int);