diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt index 375a1c097b..bb7a1452bc 100644 --- a/doc/src/Manual.txt +++ b/doc/src/Manual.txt @@ -1,7 +1,7 @@ LAMMPS Users Manual - + @@ -21,7 +21,7 @@

LAMMPS Documentation :c,h3 -7 Mar 2017 version :c,h4 +10 Mar 2017 version :c,h4 Version info: :h4 diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 5b776defca..37201275fe 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/balance.txt b/doc/src/balance.txt index 9aeb03392f..79728d6569 100644 --- a/doc/src/balance.txt +++ b/doc/src/balance.txt @@ -286,24 +286,32 @@ above. It performs a recursive coordinate bisectioning (RCB) of the simulation domain. The basic idea is as follows. The simulation domain is cut into 2 boxes by an axis-aligned cut in -the longest dimension, leaving one new box on either side of the cut. -All the processors are also partitioned into 2 groups, half assigned -to the box on the lower side of the cut, and half to the box on the -upper side. (If the processor count is odd, one side gets an extra -processor.) The cut is positioned so that the number of particles in -the lower box is exactly the number that the processors assigned to -that box should own for load balance to be perfect. This also makes -load balance for the upper box perfect. The positioning is done -iteratively, by a bisectioning method. Note that counting particles -on either side of the cut requires communication between all -processors at each iteration. +one of the dimensions, leaving one new sub-box on either side of the +cut. Which dimension is chosen for the cut depends on the particle +(weight) distribution within the parent box. Normally the longest +dimension of the box is cut, but if all (or most) of the particles are +at one end of the box, a cut may be performed in another dimension to +induce sub-boxes that are more cube-ish (3d) or square-ish (2d) in +shape. + +After the cut is made, all the processors are also partitioned into 2 +groups, half assigned to the box on the lower side of the cut, and +half to the box on the upper side. (If the processor count is odd, +one side gets an extra processor.) The cut is positioned so that the +number of (weighted) particles in the lower box is exactly the number +that the processors assigned to that box should own for load balance +to be perfect. This also makes load balance for the upper box +perfect. The positioning of the cut is done iteratively, by a +bisectioning method (median search). Note that counting particles on +either side of the cut requires communication between all processors +at each iteration. That is the procedure for the first cut. Subsequent cuts are made recursively, in exactly the same manner. The subset of processors -assigned to each box make a new cut in the longest dimension of that -box, splitting the box, the subset of processors, and the particles -in the box in two. The recursion continues until every processor is -assigned a sub-box of the entire simulation domain, and owns the +assigned to each box make a new cut in one dimension of that box, +splitting the box, the subset of processors, and the particles in the +box in two. The recursion continues until every processor is assigned +a sub-box of the entire simulation domain, and owns the (weighted) particles in that sub-box. :line diff --git a/doc/src/change_box.txt b/doc/src/change_box.txt index a41463baaf..2c7a890d4c 100644 --- a/doc/src/change_box.txt +++ b/doc/src/change_box.txt @@ -101,11 +101,11 @@ Instead you could do something like this, assuming the simulation box is non-periodic and atoms extend from 0 to 20 in all dimensions: change_box all x final -10 20 -create_atoms 1 single -5 5 5 # this will fail to insert an atom :pre +create_atoms 1 single -5 5 5 # this will fail to insert an atom :pre change_box all x final -10 20 boundary f s s create_atoms 1 single -5 5 5 -change_box boundary s s s # this will work :pre +change_box all boundary s s s # this will work :pre NOTE: Unlike the earlier "displace_box" version of this command, atom remapping is NOT performed by default. This command allows remapping diff --git a/doc/src/create_atoms.txt b/doc/src/create_atoms.txt index da9c8809d0..98c3c24a0b 100644 --- a/doc/src/create_atoms.txt +++ b/doc/src/create_atoms.txt @@ -134,6 +134,17 @@ not overlap existing atoms inappropriately, especially if molecules are being added. The "delete_atoms"_delete_atoms.html command can be used to remove overlapping atoms or molecules. +NOTE: You cannot use any of the styles explained above to create atoms +that are outside the simulation box; they will just be ignored by +LAMMPS. This is true even if you are using shrink-wrapped box +boundaries, as specified by the "boundary"_boundary.html command. +However, you can first use the "change_box"_change_box.html command to +temporarily expand the box, then add atoms via create_atoms, then +finally use change_box command again if needed to re-shrink-wrap the +new atoms. See the "change_box"_change_box.html doc page for an +example of how to do this, using the create_atoms {single} style to +insert a new atom outside the current simulation box. + :line Individual atoms are inserted by this command, unless the {mol} diff --git a/doc/src/fix_halt.txt b/doc/src/fix_halt.txt index c9295eca69..3f7650466f 100644 --- a/doc/src/fix_halt.txt +++ b/doc/src/fix_halt.txt @@ -15,15 +15,16 @@ fix ID group-ID halt N attribute operator avalue keyword value ... :pre ID, group-ID are documented in "fix"_fix.html command :ulb,l halt = style name of this fix command :l N = check halt condition every N steps :l -attribute = hstyle or v_name :l - hstyle = {bondmax} +attribute = {bondmax} or {tlimit} or v_name :l + bondmax = length of longest bond in the system + tlimit = elapsed CPU time v_name = name of "equal-style variable"_variable.html :pre operator = "<" or "<=" or ">" or ">=" or "==" or "!=" or "|^" :l avalue = numeric value to compare attribute to :l -string = text string to print with optional variable names :l zero or more keyword/value pairs may be appended :l -keyword = {error} :l - {error} value = {hard} or {soft} or {continue} :pre +keyword = {error} or {message} :l + {error} value = {hard} or {soft} or {continue} + {message} value = {yes} or {no} :pre :ule [Examples:] @@ -40,14 +41,33 @@ specified by the "run"_run.html or "minimize"_minimize.html command. The specified group-ID is ignored by this fix. -The specified {attribute} can be one of the {hstyle} options listed -above, or an "equal-style variable"_variable.html referenced as -{v_name}, where "name" is the name of a variable that has been defined -previously in the input script. +The specified {attribute} can be one of the options listed above, +namely {bondmax} or {tlimit}, or an "equal-style +variable"_variable.html referenced as {v_name}, where "name" is the +name of a variable that has been defined previously in the input +script. -The only {hstyle} option currently implemented is {bondmax}. This -will loop over all bonds in the system, compute their current -lengths, and set {attribute} to the longest bond distance. +The {bondmax} attribute will loop over all bonds in the system, +compute their current lengths, and set {attribute} to the longest bond +distance. + +The {tlimit} attribute queries the elapsed CPU time (in seconds) since +the current run began, and sets {attribute} to that value. This is an +alternative way to limit the length of a simulation run, similar to +the "timer"_timer.html timeout command. There are two differences in +using this method versus the timer command option. The first is that +the clock starts at the beginning of the current run (not when the +timer or fix command is specified), so that any setup time for the run +is not included in the elapsed time. The second is that the timer +invocation and syncing across all processors (via MPI_Allreduce) is +not performed once every {N} steps by this command. Instead it is +performed (typically) only a small number of times and the elapsed +times are used to predict when the end-of-the-run will be. Both of +these attributes can be useful when performing benchmark calculations +for a desired length of time with minmimal overhead. For example, if +a run is performing 1000s of timesteps/sec, the overhead for syncing +the timer frequently across a large number of processors may be +non-negligble. Equal-style variables evaluate to a numeric value. See the "variable"_variable.html command for a description. They calculate @@ -100,6 +120,14 @@ Note that you may wish use the "unfix"_unfix.html command on the fix halt ID, so that the same condition is not immediately triggered in a subsequent run. +The optional {message} keyword determines whether a message is printed +to the screen and logfile when the half condition is triggered. If +{message} is set to yes, a one line message with the values that +triggered the halt is printed. If {message} is set to no, no message +is printed; the run simply exits. The latter may be desirable for +post-processing tools that extract thermodyanmic information from log +files. + [Restart, fix_modify, output, run start/stop, minimize info:] No information about this fix is written to "binary restart @@ -118,4 +146,4 @@ This fix is not invoked during "energy minimization"_minimize.html. [Default:] -The option defaults are error = hard. +The option defaults are error = hard and message = yes. diff --git a/doc/src/quit.txt b/doc/src/quit.txt index 809fb2e38b..843d3de7f3 100644 --- a/doc/src/quit.txt +++ b/doc/src/quit.txt @@ -17,7 +17,7 @@ status = numerical exit status (optional) [Examples:] quit -if "$n > 10000" then "quit 1":pre +if "$n > 10000" then "quit 1" :pre [Description:] diff --git a/doc/src/thermo_style.txt b/doc/src/thermo_style.txt index e30e7023e4..36ec7bf12e 100644 --- a/doc/src/thermo_style.txt +++ b/doc/src/thermo_style.txt @@ -36,7 +36,7 @@ args = list of arguments for a particular style :l elaplong = timesteps since start of initial run in a series of runs dt = timestep size time = simulation time - cpu = elapsed CPU time in seconds + cpu = elapsed CPU time in seconds since start of this run tpcpu = time per CPU second spcpu = timesteps per CPU second cpuremain = estimated CPU time remaining in run diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index d60a9e90f7..e8c7e88324 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -1,5 +1,5 @@ -// -*- c++ -*- +// -*- c++ -*- // 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 @@ -7,6 +7,7 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" @@ -23,20 +24,31 @@ bool compare(colvar::cvc *i, colvar::cvc *j) { } -colvar::colvar(std::string const &conf) - : colvarparse(conf) +colvar::colvar() +{ + // Initialize static array once and for all + init_cv_requires(); +} + + +int colvar::init(std::string const &conf) { cvm::log("Initializing a new collective variable.\n"); + colvarparse::init(conf); + int error_code = COLVARS_OK; - get_keyval(conf, "name", this->name, - (std::string("colvar")+cvm::to_str(cvm::colvars.size()+1))); + colvarmodule *cv = cvm::main(); - if (cvm::colvar_by_name(this->name) != NULL) { + get_keyval(conf, "name", this->name, + (std::string("colvar")+cvm::to_str(cv->variables()->size()+1))); + + if ((cvm::colvar_by_name(this->name) != NULL) && + (cvm::colvar_by_name(this->name) != this)) { cvm::error("Error: this colvar cannot have the same name, \""+this->name+ "\", as another colvar.\n", INPUT_ERROR); - return; + return INPUT_ERROR; } // Initialize dependency members @@ -44,14 +56,13 @@ colvar::colvar(std::string const &conf) this->description = "colvar " + this->name; - // Initialize static array once and for all - init_cv_requires(); - kinetic_energy = 0.0; potential_energy = 0.0; error_code |= init_components(conf); - if (error_code != COLVARS_OK) return; + if (error_code != COLVARS_OK) { + return cvm::get_error(); + } size_t i; @@ -59,8 +70,6 @@ colvar::colvar(std::string const &conf) if (get_keyval(conf, "scriptedFunction", scripted_function, "", colvarparse::parse_silent)) { - // Make feature available only on user request - provide(f_cv_scripted); enable(f_cv_scripted); cvm::log("This colvar uses scripted function \"" + scripted_function + "\"."); @@ -76,8 +85,8 @@ colvar::colvar(std::string const &conf) } } if (x.type() == colvarvalue::type_notset) { - cvm::error("Could not parse scripted colvar type."); - return; + cvm::error("Could not parse scripted colvar type.", INPUT_ERROR); + return INPUT_ERROR; } cvm::log(std::string("Expecting colvar value of type ") @@ -86,8 +95,9 @@ colvar::colvar(std::string const &conf) if (x.type() == colvarvalue::type_vector) { int size; if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) { - cvm::error("Error: no size specified for vector scripted function."); - return; + cvm::error("Error: no size specified for vector scripted function.", + INPUT_ERROR); + return INPUT_ERROR; } x.vector1d_value.resize(size); } @@ -154,7 +164,7 @@ colvar::colvar(std::string const &conf) } } } - feature_states[f_cv_linear]->enabled = lin; + set_enabled(f_cv_linear, lin); } // Colvar is homogeneous if: @@ -168,7 +178,7 @@ colvar::colvar(std::string const &conf) homogeneous = false; } } - feature_states[f_cv_homogeneous]->enabled = homogeneous; + set_enabled(f_cv_homogeneous, homogeneous); } // Colvar is deemed periodic if: @@ -188,7 +198,7 @@ colvar::colvar(std::string const &conf) "Make sure that you know what you are doing!"); } } - feature_states[f_cv_periodic]->enabled = b_periodic; + set_enabled(f_cv_periodic, b_periodic); } // check that cvcs are compatible @@ -202,7 +212,7 @@ colvar::colvar(std::string const &conf) "by using components of different types. " "You must use the same type in order to " "sum them together.\n", INPUT_ERROR); - return; + return INPUT_ERROR; } } @@ -215,44 +225,110 @@ colvar::colvar(std::string const &conf) f.type(value()); f_accumulated.type(value()); + x_old.type(value()); + v_fdiff.type(value()); + v_reported.type(value()); + fj.type(value()); + ft.type(value()); + ft_reported.type(value()); + f_old.type(value()); + f_old.reset(); + + x_restart.type(value()); + after_restart = false; + reset_bias_force(); + // TODO use here information from the CVCs' own natural boundaries + error_code |= init_grid_parameters(conf); + + get_keyval(conf, "timeStepFactor", time_step_factor, 1); + + error_code |= init_extended_Lagrangian(conf); + error_code |= init_output_flags(conf); + + // Start in active state by default + enable(f_cv_active); + // Make sure dependency side-effects are correct + refresh_deps(); + + if (cvm::b_analysis) + parse_analysis(conf); + + if (cvm::debug()) + cvm::log("Done initializing collective variable \""+this->name+"\".\n"); + + return error_code; +} + + +int colvar::init_grid_parameters(std::string const &conf) +{ + colvarmodule *cv = cvm::main(); + get_keyval(conf, "width", width, 1.0); if (width <= 0.0) { cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); - return; + return INPUT_ERROR; } lower_boundary.type(value()); - lower_wall.type(value()); upper_boundary.type(value()); upper_wall.type(value()); - feature_states[f_cv_scalar]->enabled = (value().type() == colvarvalue::type_scalar); + set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar)); if (is_enabled(f_cv_scalar)) { if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) { - provide(f_cv_lower_boundary); enable(f_cv_lower_boundary); } + std::string lw_conf, uw_conf; - get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0); - if (lower_wall_k > 0.0) { + if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) { + cvm::log("Warning: lowerWallConstant and lowerWall are deprecated, " + "please define a harmonicWalls bias instead.\n"); + lower_wall.type(value()); get_keyval(conf, "lowerWall", lower_wall, lower_boundary); - enable(f_cv_lower_wall); + lw_conf = std::string("\n\ + lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\ + lowerWalls "+cvm::to_str(lower_wall)+"\n"); } if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) { - provide(f_cv_upper_boundary); enable(f_cv_upper_boundary); } - get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0); - if (upper_wall_k > 0.0) { + if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) { + cvm::log("Warning: upperWallConstant and upperWall are deprecated, " + "please define a harmonicWalls bias instead.\n"); + upper_wall.type(value()); get_keyval(conf, "upperWall", upper_wall, upper_boundary); - enable(f_cv_upper_wall); + uw_conf = std::string("\n\ + upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\ + upperWalls "+cvm::to_str(upper_wall)+"\n"); + } + + if (lw_conf.size() && uw_conf.size()) { + if (lower_wall >= upper_wall) { + cvm::error("Error: the upper wall, "+ + cvm::to_str(upper_wall)+ + ", is not higher than the lower wall, "+ + cvm::to_str(lower_wall)+".\n", + INPUT_ERROR); + return INPUT_ERROR; + } + } + + if (lw_conf.size() || uw_conf.size()) { + cvm::log("Generating a new harmonicWalls bias for compatibility purposes.\n"); + std::string const walls_conf("\n\ +harmonicWalls {\n\ + name "+this->name+"w\n\ + colvars "+this->name+"\n"+lw_conf+uw_conf+ + "}\n"); + cv->append_new_config(walls_conf); } } @@ -271,16 +347,7 @@ colvar::colvar(std::string const &conf) ", is not higher than the lower boundary, "+ cvm::to_str(lower_boundary)+".\n", INPUT_ERROR); - } - } - - if (is_enabled(f_cv_lower_wall) && is_enabled(f_cv_upper_wall)) { - if (lower_wall >= upper_wall) { - cvm::error("Error: the upper wall, "+ - cvm::to_str(upper_wall)+ - ", is not higher than the lower wall, "+ - cvm::to_str(lower_wall)+".\n", - INPUT_ERROR); + return INPUT_ERROR; } } @@ -289,83 +356,90 @@ colvar::colvar(std::string const &conf) cvm::error("Error: trying to expand boundaries that already " "cover a whole period of a periodic colvar.\n", INPUT_ERROR); + return INPUT_ERROR; } if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) { cvm::error("Error: inconsistent configuration " "(trying to expand boundaries with both " "hardLowerBoundary and hardUpperBoundary enabled).\n", INPUT_ERROR); + return INPUT_ERROR; } - get_keyval(conf, "timeStepFactor", time_step_factor, 1); + return COLVARS_OK; +} - { - bool b_extended_Lagrangian; - get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false); - if (b_extended_Lagrangian) { - cvm::real temp, tolerance, period; +int colvar::init_extended_Lagrangian(std::string const &conf) +{ + bool b_extended_Lagrangian; + get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false); - cvm::log("Enabling the extended Lagrangian term for colvar \""+ - this->name+"\".\n"); + if (b_extended_Lagrangian) { + cvm::real temp, tolerance, period; - // Make feature available only on user request - provide(f_cv_extended_Lagrangian); - enable(f_cv_extended_Lagrangian); - provide(f_cv_Langevin); + cvm::log("Enabling the extended Lagrangian term for colvar \""+ + this->name+"\".\n"); - // The extended mass will apply forces - enable(f_cv_gradient); + enable(f_cv_extended_Lagrangian); - xr.type(value()); - vr.type(value()); - fr.type(value()); + xr.type(value()); + vr.type(value()); + fr.type(value()); - const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); - if (temp <= 0.0) { - if (found) - cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR); - else - cvm::error("Error: a positive temperature must be provided, either " - "by enabling a thermostat, or through \"extendedTemp\".\n", - INPUT_ERROR); + const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); + if (temp <= 0.0) { + if (found) + cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR); + else + cvm::error("Error: a positive temperature must be provided, either " + "by enabling a thermostat, or through \"extendedTemp\".\n", + INPUT_ERROR); + return INPUT_ERROR; + } + + get_keyval(conf, "extendedFluctuation", tolerance); + if (tolerance <= 0.0) { + cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR); + return INPUT_ERROR; + } + ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); + cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2"); + + get_keyval(conf, "extendedTimeConstant", period, 200.0); + if (period <= 0.0) { + cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR); + } + ext_mass = (cvm::boltzmann() * temp * period * period) + / (4.0 * PI * PI * tolerance * tolerance); + cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)"); + + { + bool b_output_energy; + get_keyval(conf, "outputEnergy", b_output_energy, false); + if (b_output_energy) { + enable(f_cv_output_energy); } + } - get_keyval(conf, "extendedFluctuation", tolerance); - if (tolerance <= 0.0) { - cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR); - } - ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); - cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2"); - - get_keyval(conf, "extendedTimeConstant", period, 200.0); - if (period <= 0.0) { - cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR); - } - ext_mass = (cvm::boltzmann() * temp * period * period) - / (4.0 * PI * PI * tolerance * tolerance); - cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)"); - - { - bool b_output_energy; - get_keyval(conf, "outputEnergy", b_output_energy, false); - if (b_output_energy) { - enable(f_cv_output_energy); - } - } - - get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0); - if (ext_gamma < 0.0) { - cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR); - } - if (ext_gamma != 0.0) { - enable(f_cv_Langevin); - ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1 - ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt()); - } + get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0); + if (ext_gamma < 0.0) { + cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR); + return INPUT_ERROR; + } + if (ext_gamma != 0.0) { + enable(f_cv_Langevin); + ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1 + ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt()); } } + return COLVARS_OK; +} + + +int colvar::init_output_flags(std::string const &conf) +{ { bool b_output_value; get_keyval(conf, "outputValue", b_output_value, true); @@ -387,7 +461,7 @@ colvar::colvar(std::string const &conf) if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) { cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n" "The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR); - return; + return INPUT_ERROR; } } @@ -395,28 +469,7 @@ colvar::colvar(std::string const &conf) get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false); get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false); - // Start in active state by default - enable(f_cv_active); - // Make sure dependency side-effects are correct - refresh_deps(); - - x_old.type(value()); - v_fdiff.type(value()); - v_reported.type(value()); - fj.type(value()); - ft.type(value()); - ft_reported.type(value()); - f_old.type(value()); - f_old.reset(); - - x_restart.type(value()); - after_restart = false; - - if (cvm::b_analysis) - parse_analysis(conf); - - if (cvm::debug()) - cvm::log("Done initializing collective variable \""+this->name+"\".\n"); + return COLVARS_OK; } @@ -637,7 +690,7 @@ int colvar::parse_analysis(std::string const &conf) std::string runave_outfile; get_keyval(conf, "runAveOutputFile", runave_outfile, - std::string(cvm::output_prefix+"."+ + std::string(cvm::output_prefix()+"."+ this->name+".runave.traj")); size_t const this_cv_width = x.output_width(cvm::cv_width); @@ -693,7 +746,7 @@ int colvar::parse_analysis(std::string const &conf) get_keyval(conf, "corrFuncNormalize", acf_normalize, true); get_keyval(conf, "corrFuncOutputFile", acf_outfile, - std::string(cvm::output_prefix+"."+this->name+ + std::string(cvm::output_prefix()+"."+this->name+ ".corrfunc.dat")); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -730,11 +783,12 @@ colvar::~colvar() } // remove reference to this colvar from the CVM - for (std::vector::iterator cvi = cvm::colvars.begin(); - cvi != cvm::colvars.end(); + colvarmodule *cv = cvm::main(); + for (std::vector::iterator cvi = cv->variables()->begin(); + cvi != cv->variables()->end(); ++cvi) { if ( *cvi == this) { - cvm::colvars.erase(cvi); + cv->variables()->erase(cvi); break; } } @@ -892,7 +946,6 @@ int colvar::collect_cvc_values() cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n"); if (after_restart) { - after_restart = false; if (cvm::proxy->simulation_running()) { cvm::real const jump2 = dist2(x, x_restart) / (width*width); if (jump2 > 0.25) { @@ -1122,8 +1175,7 @@ int colvar::calc_colvar_properties() // initialize the restraint center in the first step to the value // just calculated from the cvcs - // TODO: put it in the restart information - if (cvm::step_relative() == 0) { + if (cvm::step_relative() == 0 && !after_restart) { xr = x; vr.reset(); // (already 0; added for clarity) } @@ -1148,6 +1200,8 @@ int colvar::calc_colvar_properties() ft_reported = ft; } + // At the end of the first update after a restart, we can reset the flag + after_restart = false; return COLVARS_OK; } @@ -1173,51 +1227,17 @@ cvm::real colvar::update_forces_energy() f -= fj; } - // Wall force - colvarvalue fw(x); - fw.reset(); - - if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) { - - if (cvm::debug()) - cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n"); - - // For a periodic colvar, both walls may be applicable at the same time - // in which case we pick the closer one - if ( (!is_enabled(f_cv_upper_wall)) || - (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) { - - cvm::real const grad = this->dist2_lgrad(x, lower_wall); - if (grad < 0.0) { - fw = -0.5 * lower_wall_k * grad; - if (cvm::debug()) - cvm::log("Applying a lower wall force("+ - cvm::to_str(fw)+") to \""+this->name+"\".\n"); - } - - } else { - - cvm::real const grad = this->dist2_lgrad(x, upper_wall); - if (grad > 0.0) { - fw = -0.5 * upper_wall_k * grad; - if (cvm::debug()) - cvm::log("Applying an upper wall force("+ - cvm::to_str(fw)+") to \""+this->name+"\".\n"); - } - } - } - // At this point f is the force f from external biases that will be applied to the // extended variable if there is one if (is_enabled(f_cv_extended_Lagrangian)) { if (cvm::debug()) { - cvm::log("Updating extended-Lagrangian degrees of freedom.\n"); + cvm::log("Updating extended-Lagrangian degree of freedom.\n"); } cvm::real dt = cvm::dt(); - colvarvalue f_ext(fr.type()); + colvarvalue f_ext(fr.type()); // force acting on the extended variable f_ext.reset(); // the total force is applied to the fictitious mass, while the @@ -1226,7 +1246,7 @@ cvm::real colvar::update_forces_energy() // 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 + wall force + // ie. spring force (fb_actual will be added just below) fr = f; f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x); @@ -1259,8 +1279,6 @@ cvm::real colvar::update_forces_energy() if (this->is_enabled(f_cv_periodic)) this->wrap(xr); } - // TODO remove the wall force - f += fw; // Now adding the force on the actual colvar (for those biases who // bypass the extended Lagrangian mass) f += fb_actual; @@ -1286,8 +1304,10 @@ cvm::real colvar::update_forces_energy() void colvar::communicate_forces() { size_t i; - if (cvm::debug()) + if (cvm::debug()) { cvm::log("Communicating forces from colvar \""+this->name+"\".\n"); + cvm::log("Force to be applied: " + cvm::to_str(f_accumulated) + "\n"); + } if (is_enabled(f_cv_scripted)) { std::vector > func_grads; @@ -1364,7 +1384,7 @@ int colvar::update_cvc_flags() active_cvc_square_norm = 0.; for (size_t i = 0; i < cvcs.size(); i++) { - cvcs[i]->feature_states[f_cvc_active]->enabled = cvc_flags[i]; + cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]); if (cvcs[i]->is_enabled()) { n_active_cvcs++; active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff; diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 2cf3d2dac5..0cbda450b8 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -227,11 +227,23 @@ public: /// Constructor - colvar(std::string const &conf); + colvar(); + + /// Main init function + int init(std::string const &conf); /// Parse the CVC configuration and allocate their data int init_components(std::string const &conf); + /// Init defaults for grid options + int init_grid_parameters(std::string const &conf); + + /// Init extended Lagrangian parameters + int init_extended_Lagrangian(std::string const &conf); + + /// Init output flags + int init_output_flags(std::string const &conf); + private: /// Parse the CVC configuration for all components of a certain type template int init_components_type(std::string const &conf, diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 48c16e887a..32cfadf3b6 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -98,7 +98,7 @@ cvm::atom_group::atom_group() cvm::atom_group::~atom_group() { - if (is_enabled(f_ag_scalable)) { + if (is_enabled(f_ag_scalable) && !b_dummy) { (cvm::proxy)->clear_atom_group(index); index = -1; } @@ -418,7 +418,7 @@ int cvm::atom_group::parse(std::string const &conf) // We need to know the fitting options to decide whether the group is scalable parse_error |= parse_fitting_options(group_conf); - if (is_available(f_ag_scalable_com) && !b_rotate) { + if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) { enable(f_ag_scalable_com); enable(f_ag_scalable); } @@ -500,14 +500,16 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf) int cvm::atom_group::add_index_group(std::string const &index_group_name) { - std::list::iterator names_i = cvm::index_group_names.begin(); - std::list >::iterator index_groups_i = cvm::index_groups.begin(); - for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) { + colvarmodule *cv = cvm::main(); + + std::list::iterator names_i = cv->index_group_names.begin(); + std::list >::iterator index_groups_i = cv->index_groups.begin(); + for ( ; names_i != cv->index_group_names.end() ; ++names_i, ++index_groups_i) { if (*names_i == index_group_name) break; } - if (names_i == cvm::index_group_names.end()) { + if (names_i == cv->index_group_names.end()) { cvm::error("Error: could not find index group "+ index_group_name+" among those provided by the index file.\n", INPUT_ERROR); diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index fdd2b6254c..3779c82aa3 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -19,20 +19,6 @@ colvarbias::colvarbias(char const *key) rank = 1; - if (bias_type == std::string("abf")) { - rank = cvm::n_abf_biases+1; - } - if (bias_type == std::string("harmonic") || - bias_type == std::string("linear")) { - rank = cvm::n_rest_biases+1; - } - if (bias_type == std::string("histogram")) { - rank = cvm::n_histo_biases+1; - } - if (bias_type == std::string("metadynamics")) { - rank = cvm::n_meta_biases+1; - } - has_data = false; b_output_energy = false; reset(); @@ -48,7 +34,11 @@ int colvarbias::init(std::string const &conf) colvarparse::init(conf); if (name.size() == 0) { + + // first initialization + cvm::log("Initializing a new \""+bias_type+"\" instance.\n"); + rank = cvm::num_biases_type(bias_type); get_keyval(conf, "name", name, bias_type+cvm::to_str(rank)); { @@ -69,7 +59,7 @@ int colvarbias::init(std::string const &conf) // lookup the associated colvars std::vector colvar_names; if (get_keyval(conf, "colvars", colvar_names)) { - if (colvars.size()) { + if (num_variables()) { cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n", INPUT_ERROR); return INPUT_ERROR; @@ -80,7 +70,7 @@ int colvarbias::init(std::string const &conf) } } - if (!colvars.size()) { + if (!num_variables()) { cvm::error("Error: no collective variables specified.\n", INPUT_ERROR); return INPUT_ERROR; } @@ -89,6 +79,8 @@ int colvarbias::init(std::string const &conf) cvm::log("Reinitializing bias \""+name+"\".\n"); } + output_prefix = cvm::output_prefix(); + get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy); return COLVARS_OK; @@ -98,7 +90,7 @@ int colvarbias::init(std::string const &conf) int colvarbias::reset() { bias_energy = 0.0; - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { colvar_forces[i].reset(); } return COLVARS_OK; @@ -132,12 +124,13 @@ int colvarbias::clear() } } + colvarmodule *cv = cvm::main(); // ...and from the colvars module - for (std::vector::iterator bi = cvm::biases.begin(); - bi != cvm::biases.end(); + for (std::vector::iterator bi = cv->biases.begin(); + bi != cv->biases.end(); ++bi) { if ( *bi == this) { - cvm::biases.erase(bi); + cv->biases.erase(bi); break; } } @@ -185,21 +178,29 @@ int colvarbias::add_colvar(std::string const &cv_name) int colvarbias::update() { - // Note: if anything is added here, it should be added also in the SMP block of calc_biases() - // TODO move here debug msg of bias update + if (cvm::debug()) { + cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n"); + } + has_data = true; + + bias_energy = 0.0; + for (size_t ir = 0; ir < num_variables(); ir++) { + colvar_forces[ir].reset(); + } + return COLVARS_OK; } void colvarbias::communicate_forces() { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { if (cvm::debug()) { cvm::log("Communicating a force to colvar \""+ - colvars[i]->name+"\".\n"); + variables(i)->name+"\".\n"); } - colvars[i]->add_bias_force(colvar_forces[i]); + variables(i)->add_bias_force(colvar_forces[i]); } } diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 12397dcb8f..6d5776d3db 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -32,17 +32,39 @@ public: /// Add a new collective variable to this bias int add_colvar(std::string const &cv_name); - /// Add a new collective variable to this bias - size_t number_of_colvars() const + /// How many variables are defined for this bias + inline size_t num_variables() const { return colvars.size(); } + /// Access the variables vector + inline std::vector *variables() + { + return &colvars; + } + + /// Access the i-th variable + inline colvar * variables(int i) const + { + return colvars[i]; + } + /// Retrieve colvar values and calculate their biasing forces /// Return bias energy virtual int update(); - // TODO: move update_bias here (share with metadynamics) + /// \brief Compute the energy of the bias with alternative values of the + /// collective variables (suitable for bias exchange) + virtual int calc_energy(std::vector const &values = + std::vector(0)) + { + cvm::error("Error: calc_energy() not implemented.\n", COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } + + /// Send forces to the collective variables + virtual void communicate_forces(); /// Load new configuration - force constant and/or centers only virtual int change_configuration(std::string const &conf); @@ -51,10 +73,13 @@ public: virtual cvm::real energy_difference(std::string const &conf); /// Give the total number of bins for a given bias. + // FIXME this is currently 1D only virtual int bin_num(); /// Calculate the bin index for a given bias. + // FIXME this is currently 1D only virtual int current_bin(); //// Give the count at a given bin index. + // FIXME this is currently 1D only virtual int bin_count(int bin_index); //// Share information between replicas, whatever it may be. virtual int replica_share(); @@ -62,9 +87,6 @@ public: /// Perform analysis tasks virtual void analyze() {} - /// Send forces to the collective variables - virtual void communicate_forces(); - /// \brief Constructor colvarbias(char const *key); @@ -135,6 +157,9 @@ public: return COLVARS_OK; } + /// Use this prefix for all output files + std::string output_prefix; + /// If this bias is communicating with other replicas through files, send it to them virtual int write_state_to_replicas() { @@ -162,7 +187,7 @@ protected: /// through each colvar object std::vector colvars; - /// \brief Current forces from this bias to the colvars + /// \brief Current forces from this bias to the variables std::vector colvar_forces; /// \brief Current energy of this bias (colvar_forces should be obtained by deriving this) diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index ccfa401697..d039004f09 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -30,9 +30,8 @@ int colvarbias_abf::init(std::string const &conf) { colvarbias::init(conf); - provide(f_cvb_scalar_variables); enable(f_cvb_scalar_variables); - provide(f_cvb_history_dependent); + enable(f_cvb_calc_pmf); // TODO relax this in case of VMD plugin if (cvm::temperature() == 0.0) @@ -221,9 +220,6 @@ colvarbias_abf::~colvarbias_abf() delete [] system_force; system_force = NULL; } - - if (cvm::n_abf_biases > 0) - cvm::n_abf_biases -= 1; } @@ -319,11 +315,11 @@ int colvarbias_abf::update() } // update the output prefix; TODO: move later to setup_output() function - if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) { - // This is the only ABF bias - output_prefix = cvm::output_prefix; + if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) { + // This is the only bias computing PMFs + output_prefix = cvm::output_prefix(); } else { - output_prefix = cvm::output_prefix + "." + this->name; + output_prefix = cvm::output_prefix() + "." + this->name; } if (output_freq && (cvm::step_absolute() % output_freq) == 0) { diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 9e414d4db4..d096ac3daf 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -40,11 +40,8 @@ int colvarbias_alb::init(std::string const &conf) { colvarbias::init(conf); - provide(f_cvb_scalar_variables); enable(f_cvb_scalar_variables); - provide(f_cvb_history_dependent); - size_t i; // get the initial restraint centers @@ -138,8 +135,6 @@ int colvarbias_alb::init(std::string const &conf) colvarbias_alb::~colvarbias_alb() { - if (cvm::n_rest_biases > 0) - cvm::n_rest_biases -= 1; } diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index 900ad213d5..502a7455b1 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -24,10 +24,7 @@ int colvarbias_histogram::init(std::string const &conf) { colvarbias::init(conf); - provide(f_cvb_scalar_variables); enable(f_cvb_scalar_variables); - - provide(f_cvb_history_dependent); enable(f_cvb_history_dependent); size_t i; @@ -104,9 +101,6 @@ colvarbias_histogram::~colvarbias_histogram() delete grid; grid = NULL; } - - if (cvm::n_histo_biases > 0) - cvm::n_histo_biases -= 1; } @@ -127,14 +121,14 @@ int colvarbias_histogram::update() // At the first timestep, we need to assign out_name since // output_prefix is unset during the constructor if (cvm::step_relative() == 0) { - out_name = cvm::output_prefix + "." + this->name + ".dat"; + out_name = cvm::output_prefix() + "." + this->name + ".dat"; cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\""); } } if (out_name_dx.size() == 0) { if (cvm::step_relative() == 0) { - out_name_dx = cvm::output_prefix + "." + this->name + ".dx"; + out_name_dx = cvm::output_prefix() + "." + this->name + ".dx"; cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\""); } } diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 01a04d1a01..b0acfe974a 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -43,7 +43,7 @@ int colvarbias_meta::init(std::string const &conf) { colvarbias::init(conf); - provide(f_cvb_history_dependent); + enable(f_cvb_calc_pmf); get_keyval(conf, "hillWeight", hill_weight, 0.0); if (hill_weight > 0.0) { @@ -59,9 +59,9 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0); cvm::log("Half-widths of the Gaussian hills (sigma's):\n"); - for (size_t i = 0; i < colvars.size(); i++) { - cvm::log(colvars[i]->name+std::string(": ")+ - cvm::to_str(0.5 * colvars[i]->width * hill_width)); + for (size_t i = 0; i < num_variables(); i++) { + cvm::log(variables(i)->name+std::string(": ")+ + cvm::to_str(0.5 * variables(i)->width * hill_width)); } { @@ -73,8 +73,10 @@ int colvarbias_meta::init(std::string const &conf) comm = single_replica; } - // This implies gradients for all colvars - enable(f_cvb_apply_force); + // in all cases, the first replica is this bias itself + if (replicas.size() == 0) { + replicas.push_back(this); + } get_keyval(conf, "useGrids", use_grids, true); @@ -84,14 +86,14 @@ int colvarbias_meta::init(std::string const &conf) expand_grids = false; size_t i; - for (i = 0; i < colvars.size(); i++) { - colvars[i]->enable(f_cv_grid); - if (colvars[i]->expand_boundaries) { + for (i = 0; i < num_variables(); i++) { + variables(i)->enable(f_cv_grid); + if (variables(i)->expand_boundaries) { expand_grids = true; cvm::log("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": Will expand grids when the colvar \""+ - colvars[i]->name+"\" approaches its boundaries.\n"); + variables(i)->name+"\" approaches its boundaries.\n"); } } @@ -100,7 +102,7 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent); if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) { cvm::log("Option \"saveFreeEnergyFile\" is deprecated, " - "please use \"keepFreeEnergyFiles\" instead."); + "please use \"keepFreeEnergyFile\" instead."); } get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save); @@ -154,6 +156,20 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false); + init_well_tempered_params(conf); + init_ebmeta_params(conf); + + if (cvm::debug()) + cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); + + save_delimiters = false; + return COLVARS_OK; +} + + +int colvarbias_meta::init_well_tempered_params(std::string const &conf) +{ // for well-tempered metadynamics get_keyval(conf, "wellTempered", well_tempered, false); get_keyval(conf, "biasTemperature", bias_temperature, -1.0); @@ -164,8 +180,12 @@ int colvarbias_meta::init(std::string const &conf) cvm::log("Well-tempered metadynamics is used.\n"); cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n"); } + return COLVARS_OK; +} +int colvarbias_meta::init_ebmeta_params(std::string const &conf) +{ // for ebmeta target_dist = NULL; get_keyval(conf, "ebMeta", ebmeta, false); @@ -203,11 +223,6 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0); } - if (cvm::debug()) - cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); - - save_delimiters = false; return COLVARS_OK; } @@ -234,9 +249,6 @@ colvarbias_meta::~colvarbias_meta() delete target_dist; target_dist = NULL; } - - if (cvm::n_meta_biases > 0) - cvm::n_meta_biases -= 1; } @@ -314,23 +326,45 @@ colvarbias_meta::delete_hill(hill_iter &h) int colvarbias_meta::update() { - if (cvm::debug()) - cvm::log("Updating the metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); + int error_code = COLVARS_OK; + // update base class + error_code |= colvarbias::update(); + + // update grid definition, if needed + error_code |= update_grid_params(); + // add new biasing energy/forces + error_code |= update_bias(); + // update grid content to reflect new bias + error_code |= update_grid_data(); + + if (comm != single_replica && + (cvm::step_absolute() % replica_update_freq) == 0) { + // sync with the other replicas (if needed) + error_code |= replica_share(); + } + + error_code |= calc_energy(); + error_code |= calc_forces(); + + return error_code; +} + + +int colvarbias_meta::update_grid_params() +{ if (use_grids) { std::vector curr_bin = hills_energy->get_colvars_index(); + if (cvm::debug()) { + cvm::log("Metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": current coordinates on the grid: "+ + cvm::to_str(curr_bin)+".\n"); + } if (expand_grids) { - // first of all, expand the grids, if specified - if (cvm::debug()) - cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": current coordinates on the grid: "+ - cvm::to_str(curr_bin)+".\n"); - bool changed_grids = false; int const min_buffer = (3 * (size_t) std::floor(hill_width)) + 1; @@ -339,9 +373,9 @@ int colvarbias_meta::update() std::vector new_lower_boundaries(hills_energy->lower_boundaries); std::vector new_upper_boundaries(hills_energy->upper_boundaries); - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { - if (! colvars[i]->expand_boundaries) + if (! variables(i)->expand_boundaries) continue; cvm::real &new_lb = new_lower_boundaries[i].real_value; @@ -349,10 +383,10 @@ int colvarbias_meta::update() int &new_size = new_sizes[i]; bool changed_lb = false, changed_ub = false; - if (!colvars[i]->hard_lower_boundary) + if (!variables(i)->hard_lower_boundary) if (curr_bin[i] < min_buffer) { int const extra_points = (min_buffer - curr_bin[i]); - new_lb -= extra_points * colvars[i]->width; + new_lb -= extra_points * variables(i)->width; new_size += extra_points; // changed offset in this direction => the pointer needs to // be changed, too @@ -362,21 +396,21 @@ int colvarbias_meta::update() cvm::log("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": new lower boundary for colvar \""+ - colvars[i]->name+"\", at "+ + variables(i)->name+"\", at "+ cvm::to_str(new_lower_boundaries[i])+".\n"); } - if (!colvars[i]->hard_upper_boundary) + if (!variables(i)->hard_upper_boundary) if (curr_bin[i] > new_size - min_buffer - 1) { int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer); - new_ub += extra_points * colvars[i]->width; + new_ub += extra_points * variables(i)->width; new_size += extra_points; changed_ub = true; cvm::log("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": new upper boundary for colvar \""+ - colvars[i]->name+"\", at "+ + variables(i)->name+"\", at "+ cvm::to_str(new_upper_boundaries[i])+".\n"); } @@ -401,7 +435,7 @@ int colvarbias_meta::update() new_hills_energy_gradients->lower_boundaries = new_lower_boundaries; new_hills_energy_gradients->upper_boundaries = new_upper_boundaries; - new_hills_energy_gradients->setup(new_sizes, 0.0, colvars.size()); + new_hills_energy_gradients->setup(new_sizes, 0.0, num_variables()); new_hills_energy->map_grid(*hills_energy); new_hills_energy_gradients->map_grid(*hills_energy_gradients); @@ -418,25 +452,32 @@ int colvarbias_meta::update() } } } + return COLVARS_OK; +} + +int colvarbias_meta::update_bias() +{ // add a new hill if the required time interval has passed - if ((cvm::step_absolute() % new_hill_freq) == 0) { + if ((cvm::step_absolute() % new_hill_freq) == 0 && + is_enabled(f_cvb_history_dependent)) { - if (cvm::debug()) + if (cvm::debug()) { cvm::log("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n"); + } cvm::real hills_scale=1.0; if (ebmeta) { - hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index()); - if(cvm::step_absolute() <= long(ebmeta_equil_steps)) { - cvm::real const hills_lambda = - (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) / - (cvm::real(ebmeta_equil_steps)); - hills_scale = hills_lambda + (1-hills_lambda)*hills_scale; - } + hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index()); + if(cvm::step_absolute() <= long(ebmeta_equil_steps)) { + cvm::real const hills_lambda = + (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) / + (cvm::real(ebmeta_equil_steps)); + hills_scale = hills_lambda + (1-hills_lambda)*hills_scale; + } } if (well_tempered) { @@ -471,160 +512,165 @@ int colvarbias_meta::update() } } - // sync with the other replicas (if needed) - if (comm != single_replica) { + return COLVARS_OK; +} - // reread the replicas registry - if ((cvm::step_absolute() % replica_update_freq) == 0) { - update_replicas_registry(); - // empty the output buffer - if (replica_hills_os.is_open()) - replica_hills_os.flush(); - read_replica_files(); +int colvarbias_meta::update_grid_data() +{ + if ((cvm::step_absolute() % grids_freq) == 0) { + // map the most recent gaussians to the grids + project_hills(new_hills_begin, hills.end(), + hills_energy, hills_energy_gradients); + new_hills_begin = hills.end(); + + // TODO: we may want to condense all into one replicas array, + // including "this" as the first element + if (comm == multiple_replicas) { + for (size_t ir = 0; ir < replicas.size(); ir++) { + replicas[ir]->project_hills(replicas[ir]->new_hills_begin, + replicas[ir]->hills.end(), + replicas[ir]->hills_energy, + replicas[ir]->hills_energy_gradients); + replicas[ir]->new_hills_begin = replicas[ir]->hills.end(); + } } } - // calculate the biasing energy and forces - bias_energy = 0.0; - for (size_t ir = 0; ir < colvars.size(); ir++) { - colvar_forces[ir].reset(); + return COLVARS_OK; +} + + +int colvarbias_meta::calc_energy(std::vector const &values) +{ + size_t ir = 0; + + for (ir = 0; ir < replicas.size(); ir++) { + replicas[ir]->bias_energy = 0.0; } - if (comm == multiple_replicas) - for (size_t ir = 0; ir < replicas.size(); ir++) { - replicas[ir]->bias_energy = 0.0; - for (size_t ic = 0; ic < colvars.size(); ic++) { - replicas[ir]->colvar_forces[ic].reset(); + + std::vector const curr_bin = values.size() ? + hills_energy->get_colvars_index(values) : + hills_energy->get_colvars_index(); + + if (hills_energy->index_ok(curr_bin)) { + // index is within the grid: get the energy from there + for (ir = 0; ir < replicas.size(); ir++) { + + bias_energy += replicas[ir]->hills_energy->value(curr_bin); + if (cvm::debug()) { + cvm::log("Metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": current coordinates on the grid: "+ + cvm::to_str(curr_bin)+".\n"); + cvm::log("Grid energy = "+cvm::to_str(bias_energy)+".\n"); } } - - if (use_grids) { - - // get the forces from the grid - - if ((cvm::step_absolute() % grids_freq) == 0) { - // map the most recent gaussians to the grids - project_hills(new_hills_begin, hills.end(), - hills_energy, hills_energy_gradients); - new_hills_begin = hills.end(); - - // TODO: we may want to condense all into one replicas array, - // including "this" as the first element - if (comm == multiple_replicas) { - for (size_t ir = 0; ir < replicas.size(); ir++) { - replicas[ir]->project_hills(replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - replicas[ir]->hills_energy, - replicas[ir]->hills_energy_gradients); - replicas[ir]->new_hills_begin = replicas[ir]->hills.end(); - } - } - } - - std::vector curr_bin = hills_energy->get_colvars_index(); - if (cvm::debug()) - cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": current coordinates on the grid: "+ - cvm::to_str(curr_bin)+".\n"); - - if (hills_energy->index_ok(curr_bin)) { - - // within the grid: add the energy and the forces from there - - bias_energy += hills_energy->value(curr_bin); - for (size_t ic = 0; ic < colvars.size(); ic++) { - cvm::real const *f = &(hills_energy_gradients->value(curr_bin)); - colvar_forces[ic].real_value += -1.0 * f[ic]; - // the gradients are stored, not the forces - } - if (comm == multiple_replicas) - for (size_t ir = 0; ir < replicas.size(); ir++) { - bias_energy += replicas[ir]->hills_energy->value(curr_bin); - cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin)); - for (size_t ic = 0; ic < colvars.size(); ic++) { - colvar_forces[ic].real_value += -1.0 * f[ic]; - } - } - - } else { - - // off the grid: compute analytically only the hills at the grid's edges - - calc_hills(hills_off_grid.begin(), hills_off_grid.end(), bias_energy); - for (size_t ic = 0; ic < colvars.size(); ic++) { - calc_hills_force(ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces); - } - - if (comm == multiple_replicas) - for (size_t ir = 0; ir < replicas.size(); ir++) { - calc_hills(replicas[ir]->hills_off_grid.begin(), - replicas[ir]->hills_off_grid.end(), - bias_energy); - for (size_t ic = 0; ic < colvars.size(); ic++) { - calc_hills_force(ic, - replicas[ir]->hills_off_grid.begin(), - replicas[ir]->hills_off_grid.end(), - colvar_forces); - } - } + } else { + // off the grid: compute analytically only the hills at the grid's edges + for (ir = 0; ir < replicas.size(); ir++) { + calc_hills(replicas[ir]->hills_off_grid.begin(), + replicas[ir]->hills_off_grid.end(), + bias_energy, + values); } } // now include the hills that have not been binned yet (starting // from new_hills_begin) - calc_hills(new_hills_begin, hills.end(), bias_energy); - for (size_t ic = 0; ic < colvars.size(); ic++) { - calc_hills_force(ic, new_hills_begin, hills.end(), colvar_forces); - } - - if (cvm::debug()) - cvm::log("Hills energy = "+cvm::to_str(bias_energy)+ - ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); - - if (cvm::debug()) - cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": adding the forces from the other replicas.\n"); - - if (comm == multiple_replicas) - for (size_t ir = 0; ir < replicas.size(); ir++) { - calc_hills(replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - bias_energy); - for (size_t ic = 0; ic < colvars.size(); ic++) { - calc_hills_force(ic, - replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - colvar_forces); - } - if (cvm::debug()) - cvm::log("Hills energy = "+cvm::to_str(bias_energy)+ - ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); + for (ir = 0; ir < replicas.size(); ir++) { + calc_hills(replicas[ir]->new_hills_begin, + replicas[ir]->hills.end(), + bias_energy); + if (cvm::debug()) { + cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n"); } + } return COLVARS_OK; } +int colvarbias_meta::calc_forces(std::vector const &values) +{ + size_t ir = 0, ic = 0; + for (ir = 0; ir < replicas.size(); ir++) { + for (ic = 0; ic < num_variables(); ic++) { + replicas[ir]->colvar_forces[ic].reset(); + } + } + + std::vector const curr_bin = values.size() ? + hills_energy->get_colvars_index(values) : + hills_energy->get_colvars_index(); + + if (hills_energy->index_ok(curr_bin)) { + for (ir = 0; ir < replicas.size(); ir++) { + cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin)); + for (ic = 0; ic < num_variables(); ic++) { + // the gradients are stored, not the forces + colvar_forces[ic].real_value += -1.0 * f[ic]; + } + } + } else { + // off the grid: compute analytically only the hills at the grid's edges + for (ir = 0; ir < replicas.size(); ir++) { + for (ic = 0; ic < num_variables(); ic++) { + calc_hills_force(ic, + replicas[ir]->hills_off_grid.begin(), + replicas[ir]->hills_off_grid.end(), + colvar_forces, + values); + } + } + } + + // now include the hills that have not been binned yet (starting + // from new_hills_begin) + + if (cvm::debug()) { + cvm::log("Metadynamics bias \""+this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": adding the forces from the other replicas.\n"); + } + + for (ir = 0; ir < replicas.size(); ir++) { + for (ic = 0; ic < num_variables(); ic++) { + calc_hills_force(ic, + replicas[ir]->new_hills_begin, + replicas[ir]->hills.end(), + colvar_forces, + values); + if (cvm::debug()) { + cvm::log("Hills forces = "+cvm::to_str(colvar_forces)+".\n"); + } + } + } + + return COLVARS_OK; +} + + + void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, cvm::real &energy, std::vector const &colvar_values) { - std::vector curr_values(colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - curr_values[i].type(colvars[i]->value()); + size_t i = 0; + std::vector curr_values(num_variables()); + for (i = 0; i < num_variables(); i++) { + curr_values[i].type(variables(i)->value()); } if (colvar_values.size()) { - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { curr_values[i] = colvar_values[i]; } } else { - for (size_t i = 0; i < colvars.size(); i++) { - curr_values[i] = colvars[i]->value(); + for (i = 0; i < num_variables(); i++) { + curr_values[i] = variables(i)->value(); } } @@ -632,11 +678,11 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, // compute the gaussian exponent cvm::real cv_sqdev = 0.0; - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { colvarvalue const &x = curr_values[i]; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; - cv_sqdev += (colvars[i]->dist2(x, center)) / (half_width*half_width); + cv_sqdev += (variables(i)->dist2(x, center)) / (half_width*half_width); } // compute the gaussian @@ -658,14 +704,14 @@ void colvarbias_meta::calc_hills_force(size_t const &i, std::vector const &values) { // Retrieve the value of the colvar - colvarvalue const x(values.size() ? values[i] : colvars[i]->value()); + colvarvalue const x(values.size() ? values[i] : variables(i)->value()); // do the type check only once (all colvarvalues in the hills series // were already saved with their types matching those in the // colvars) hill_iter h; - switch (colvars[i]->value().type()) { + switch (variables(i)->value().type()) { case colvarvalue::type_scalar: for (h = h_first; h != h_last; h++) { @@ -674,7 +720,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i, cvm::real const half_width = 0.5 * h->widths[i]; forces[i].real_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * - (colvars[i]->dist2_lgrad(x, center)).real_value ); + (variables(i)->dist2_lgrad(x, center)).real_value ); } break; @@ -687,7 +733,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i, cvm::real const half_width = 0.5 * h->widths[i]; forces[i].rvector_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * - (colvars[i]->dist2_lgrad(x, center)).rvector_value ); + (variables(i)->dist2_lgrad(x, center)).rvector_value ); } break; @@ -699,7 +745,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i, cvm::real const half_width = 0.5 * h->widths[i]; forces[i].quaternion_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * - (colvars[i]->dist2_lgrad(x, center)).quaternion_value ); + (variables(i)->dist2_lgrad(x, center)).quaternion_value ); } break; @@ -710,7 +756,7 @@ void colvarbias_meta::calc_hills_force(size_t const &i, cvm::real const half_width = 0.5 * h->widths[i]; forces[i].vector1d_value += ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * - (colvars[i]->dist2_lgrad(x, center)).vector1d_value ); + (variables(i)->dist2_lgrad(x, center)).vector1d_value ); } break; @@ -739,13 +785,13 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, // TODO: improve it by looping over a small subgrid instead of the whole grid - std::vector colvar_values(colvars.size()); - std::vector colvar_forces_scalar(colvars.size()); + std::vector colvar_values(num_variables()); + std::vector colvar_forces_scalar(num_variables()); std::vector he_ix = he->new_index(); std::vector hg_ix = (hg != NULL) ? hg->new_index() : std::vector (0); cvm::real hills_energy_here = 0.0; - std::vector hills_forces_here(colvars.size(), 0.0); + std::vector hills_forces_here(num_variables(), 0.0); size_t count = 0; size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1))); @@ -757,7 +803,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, (he->index_ok(he_ix)) && (hg->index_ok(hg_ix)); count++) { size_t i; - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i); } @@ -766,7 +812,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, calc_hills(h_first, h_last, hills_energy_here, colvar_values); he->acc_value(he_ix, hills_energy_here); - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { hills_forces_here[i].reset(); calc_hills_force(i, h_first, h_last, hills_forces_here, colvar_values); colvar_forces_scalar[i] = hills_forces_here[i].real_value; @@ -795,7 +841,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, for ( ; (he->index_ok(he_ix)); ) { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i); } @@ -851,6 +897,21 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first // ********************************************************************** +int colvarbias_meta::replica_share() +{ + // sync with the other replicas (if needed) + if (comm == multiple_replicas) { + // reread the replicas registry + update_replicas_registry(); + // empty the output buffer + if (replica_hills_os.is_open()) + replica_hills_os.flush(); + read_replica_files(); + } + return COLVARS_OK; +} + + void colvarbias_meta::update_replicas_registry() { if (cvm::debug()) @@ -975,7 +1036,6 @@ void colvarbias_meta::update_replicas_registry() } } - if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+ cvm::to_str(replicas.size())+" elements.\n"); @@ -984,7 +1044,8 @@ void colvarbias_meta::update_replicas_registry() void colvarbias_meta::read_replica_files() { - for (size_t ir = 0; ir < replicas.size(); ir++) { + // Note: we start from the 2nd replica. + for (size_t ir = 1; ir < replicas.size(); ir++) { if (! (replicas[ir])->replica_state_file_in_sync) { // if a new state file is being read, the hills file is also new @@ -1352,9 +1413,9 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) cvm::real h_weight; get_keyval(data, "weight", h_weight, hill_weight, parse_silent); - std::vector h_centers(colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - h_centers[i].type(colvars[i]->value()); + std::vector h_centers(num_variables()); + for (size_t i = 0; i < num_variables(); i++) { + h_centers[i].type(variables(i)->value()); } { // it is safer to read colvarvalue objects one at a time; @@ -1362,14 +1423,14 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) std::string centers_input; key_lookup(data, "centers", centers_input); std::istringstream centers_is(centers_input); - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { centers_is >> h_centers[i]; } } - std::vector h_widths(colvars.size()); + std::vector h_widths(num_variables()); get_keyval(data, "widths", h_widths, - std::vector (colvars.size(), (std::sqrt(2.0 * PI) / 2.0)), + std::vector(num_variables(), (std::sqrt(2.0 * PI) / 2.0)), parse_silent); std::string h_replica = ""; @@ -1406,6 +1467,13 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) int colvarbias_meta::setup_output() { + output_prefix = cvm::output_prefix(); + if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) { + // if this is not the only free energy integrator, append + // this bias's name, to distinguish it from the output of the other + // biases producing a .pmf file + output_prefix += ("."+this->name); + } if (comm == multiple_replicas) { @@ -1421,10 +1489,10 @@ int colvarbias_meta::setup_output() // those by another replica replica_hills_file = (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills"); + cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills"); replica_state_file = (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state"); + cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state"); delete[] pwd; // now register this replica @@ -1492,13 +1560,14 @@ int colvarbias_meta::setup_output() } if (b_hills_traj) { - std::string const traj_file_name(cvm::output_prefix+ + std::string const traj_file_name(cvm::output_prefix()+ ".colvars."+this->name+ ( (comm != single_replica) ? ("."+replica_id) : ("") )+ ".hills.traj"); if (!hills_traj_os.is_open()) { + cvm::backup_file(traj_file_name.c_str()); hills_traj_os.open(traj_file_name.c_str()); } if (!hills_traj_os.is_open()) @@ -1585,16 +1654,6 @@ void colvarbias_meta::write_pmf() colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy); pmf->setup(); - std::string fes_file_name_prefix(cvm::output_prefix); - - if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) { - // if this is not the only free energy integrator, append - // this bias's name, to distinguish it from the output of the other - // biases producing a .pmf file - // TODO: fix for ABF with updateBias == no - fes_file_name_prefix += ("."+this->name); - } - if ((comm == single_replica) || (dump_replica_fes)) { // output the PMF from this instance or replica pmf->reset(); @@ -1607,7 +1666,7 @@ void colvarbias_meta::write_pmf() pmf->multiply_constant(well_temper_scale); } { - std::string const fes_file_name(fes_file_name_prefix + + std::string const fes_file_name(this->output_prefix + ((comm != single_replica) ? ".partial" : "") + (dump_fes_save ? "."+cvm::to_str(cvm::step_absolute()) : "") + @@ -1632,7 +1691,7 @@ void colvarbias_meta::write_pmf() cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature; pmf->multiply_constant(well_temper_scale); } - std::string const fes_file_name(fes_file_name_prefix + + std::string const fes_file_name(this->output_prefix + (dump_fes_save ? "."+cvm::to_str(cvm::step_absolute()) : "") + ".pmf"); diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index a88a34ba00..01921eaf64 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -36,8 +36,20 @@ public: colvarbias_meta(char const *key); virtual int init(std::string const &conf); + virtual int init_well_tempered_params(std::string const &conf); + virtual int init_ebmeta_params(std::string const &conf); virtual ~colvarbias_meta(); + virtual int update(); + virtual int update_grid_params(); + virtual int update_bias(); + virtual int update_grid_data(); + virtual int replica_share(); + + virtual int calc_energy(std::vector const &values = + std::vector(0)); + virtual int calc_forces(std::vector const &values = + std::vector(0)); virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &state_conf); @@ -102,18 +114,18 @@ protected: /// \brief Calculate the values of the hills, incrementing /// bias_energy virtual void calc_hills(hill_iter h_first, - hill_iter h_last, - cvm::real &energy, - std::vector const &values = std::vector (0)); + hill_iter h_last, + cvm::real &energy, + std::vector const &values = std::vector(0)); /// \brief Calculate the forces acting on the i-th colvar, /// incrementing colvar_forces[i]; must be called after calc_hills /// each time the values of the colvars are changed virtual void calc_hills_force(size_t const &i, - hill_iter h_first, - hill_iter h_last, - std::vector &forces, - std::vector const &values = std::vector (0)); + hill_iter h_first, + hill_iter h_last, + std::vector &forces, + std::vector const &values = std::vector(0)); /// Height of new hills diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 84630984e5..159d9eae64 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -33,17 +33,15 @@ int colvarbias_restraint::init(std::string const &conf) int colvarbias_restraint::update() { - bias_energy = 0.0; - - if (cvm::debug()) - cvm::log("Updating the restraint bias \""+this->name+"\".\n"); + // Update base class (bias_energy and colvar_forces are zeroed there) + colvarbias::update(); // Force and energy calculation - for (size_t i = 0; i < colvars.size(); i++) { - colvar_forces[i].type(colvars[i]->value()); + for (size_t i = 0; i < num_variables(); i++) { + bias_energy += restraint_potential(i); + colvar_forces[i].type(variables(i)->value()); colvar_forces[i].is_derivative(); colvar_forces[i] = restraint_force(i); - bias_energy += restraint_potential(i); } if (cvm::debug()) @@ -59,8 +57,6 @@ int colvarbias_restraint::update() colvarbias_restraint::~colvarbias_restraint() { - if (cvm::n_rest_biases > 0) - cvm::n_rest_biases -= 1; } @@ -102,18 +98,18 @@ int colvarbias_restraint_centers::init(std::string const &conf) bool null_centers = (colvar_centers.size() == 0); if (null_centers) { // try to initialize the restraint centers for the first time - colvar_centers.resize(colvars.size()); - colvar_centers_raw.resize(colvars.size()); - for (i = 0; i < colvars.size(); i++) { - colvar_centers[i].type(colvars[i]->value()); + colvar_centers.resize(num_variables()); + colvar_centers_raw.resize(num_variables()); + for (i = 0; i < num_variables(); i++) { + colvar_centers[i].type(variables(i)->value()); colvar_centers[i].reset(); - colvar_centers_raw[i].type(colvars[i]->value()); + colvar_centers_raw[i].type(variables(i)->value()); colvar_centers_raw[i].reset(); } } if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { if (cvm::debug()) { cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n"); } @@ -129,7 +125,7 @@ int colvarbias_restraint_centers::init(std::string const &conf) return INPUT_ERROR; } - if (colvar_centers.size() != colvars.size()) { + if (colvar_centers.size() != num_variables()) { cvm::error("Error: number of centers does not match " "that of collective variables.\n", INPUT_ERROR); return INPUT_ERROR; @@ -142,10 +138,10 @@ int colvarbias_restraint_centers::init(std::string const &conf) int colvarbias_restraint_centers::change_configuration(std::string const &conf) { if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].type(colvars[i]->value()); + for (size_t i = 0; i < num_variables(); i++) { + colvar_centers[i].type(variables(i)->value()); colvar_centers[i].apply_constraints(); - colvar_centers_raw[i].type(colvars[i]->value()); + colvar_centers_raw[i].type(variables(i)->value()); colvar_centers_raw[i] = colvar_centers[i]; } } @@ -269,7 +265,7 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf) size_t i; if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { - if (colvar_centers.size() != colvars.size()) { + if (colvar_centers.size() != num_variables()) { cvm::error("Error: number of target centers does not match " "that of collective variables.\n"); } @@ -308,9 +304,9 @@ int colvarbias_restraint_centers_moving::update() // at each simulation step (or stage, if applicable) // (take current stage into account: it can be non-zero // if we are restarting a staged calculation) - centers_incr.resize(colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - centers_incr[i].type(colvars[i]->value()); + centers_incr.resize(num_variables()); + for (size_t i = 0; i < num_variables(); i++) { + centers_incr[i].type(variables(i)->value()); centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / cvm::real( target_nstages ? (target_nstages - stage) : (target_nsteps - cvm::step_absolute())); @@ -326,10 +322,10 @@ int colvarbias_restraint_centers_moving::update() && (cvm::step_absolute() % target_nsteps) == 0 && stage < target_nstages) { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { colvar_centers_raw[i] += centers_incr[i]; colvar_centers[i] = colvar_centers_raw[i]; - colvars[i]->wrap(colvar_centers[i]); + variables(i)->wrap(colvar_centers[i]); colvar_centers[i].apply_constraints(); } stage++; @@ -341,10 +337,10 @@ int colvarbias_restraint_centers_moving::update() } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) { // move the restraint centers in the direction of the targets // (slow growth) - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { colvar_centers_raw[i] += centers_incr[i]; colvar_centers[i] = colvar_centers_raw[i]; - colvars[i]->wrap(colvar_centers[i]); + variables(i)->wrap(colvar_centers[i]); colvar_centers[i].apply_constraints(); } } @@ -363,7 +359,7 @@ int colvarbias_restraint_centers_moving::update_acc_work() { if (b_output_acc_work) { if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { // project forces on the calculated increments at this step acc_work += colvar_forces[i] * centers_incr[i]; } @@ -381,14 +377,14 @@ std::string const colvarbias_restraint_centers_moving::get_state_params() const if (b_chg_centers) { size_t i; os << "centers "; - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << colvar_centers[i]; } os << "\n"; os << "centers_raw "; - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << colvar_centers_raw[i]; @@ -429,10 +425,10 @@ int colvarbias_restraint_centers_moving::set_state_params(std::string const &con std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os) { if (b_output_centers) { - for (size_t i = 0; i < colvars.size(); i++) { - size_t const this_cv_width = (colvars[i]->value()).output_width(cvm::cv_width); + for (size_t i = 0; i < num_variables(); i++) { + size_t const this_cv_width = (variables(i)->value()).output_width(cvm::cv_width); os << " x0_" - << cvm::wrap_string(colvars[i]->name, this_cv_width-3); + << cvm::wrap_string(variables(i)->name, this_cv_width-3); } } @@ -448,7 +444,7 @@ std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostrea std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os) { if (b_output_centers) { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << colvar_centers[i]; @@ -539,9 +535,9 @@ int colvarbias_restraint_k_moving::update() } force_k = starting_force_k + (target_force_k - starting_force_k) * std::pow(lambda, force_k_exp); - cvm::log("Restraint " + this->name + ", stage " + - cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); - cvm::log("Setting force constant to " + cvm::to_str(force_k)); + cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + + " : lambda = " + cvm::to_str(lambda) + + ", k = " + cvm::to_str(force_k)); } // TI calculation: estimate free energy derivative @@ -557,7 +553,7 @@ int colvarbias_restraint_k_moving::update() // Square distance normalized by square colvar width cvm::real dist_sq = 0.0; - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { dist_sq += d_restraint_potential_dk(i); } @@ -569,7 +565,8 @@ int colvarbias_restraint_k_moving::update() if (cvm::step_absolute() % target_nsteps == 0 && cvm::step_absolute() > 0) { - cvm::log("Lambda= " + cvm::to_str(lambda) + " dA/dLambda= " + cvm::log("Restraint " + this->name + " Lambda= " + + cvm::to_str(lambda) + " dA/dLambda= " + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps))); // ...and move on to the next one @@ -584,9 +581,9 @@ int colvarbias_restraint_k_moving::update() } force_k = starting_force_k + (target_force_k - starting_force_k) * std::pow(lambda, force_k_exp); - cvm::log("Restraint " + this->name + ", stage " + - cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); - cvm::log("Setting force constant to " + cvm::to_str(force_k)); + cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + + " : lambda = " + cvm::to_str(lambda) + + ", k = " + cvm::to_str(force_k)); } } @@ -721,11 +718,11 @@ int colvarbias_restraint_harmonic::init(std::string const &conf) colvarbias_restraint_centers_moving::init(conf); colvarbias_restraint_k_moving::init(conf); - for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->width != 1.0) - cvm::log("The force constant for colvar \""+colvars[i]->name+ + 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 / (colvars[i]->width * colvars[i]->width))+ + cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+ " according to the specified width.\n"); } @@ -751,22 +748,22 @@ int colvarbias_restraint_harmonic::update() cvm::real colvarbias_restraint_harmonic::restraint_potential(size_t i) const { - return 0.5 * force_k / (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]); + return 0.5 * force_k / (variables(i)->width * variables(i)->width) * + variables(i)->dist2(variables(i)->value(), colvar_centers[i]); } colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const { - return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2_lgrad(colvars[i]->value(), colvar_centers[i]); + return -0.5 * force_k / (variables(i)->width * variables(i)->width) * + variables(i)->dist2_lgrad(variables(i)->value(), colvar_centers[i]); } cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const { - return 0.5 / (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]); + return 0.5 / (variables(i)->width * variables(i)->width) * + variables(i)->dist2(variables(i)->value(), colvar_centers[i]); } @@ -840,6 +837,8 @@ colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char co colvarbias_restraint_moving(key), colvarbias_restraint_k_moving(key) { + lower_wall_k = 0.0; + upper_wall_k = 0.0; } @@ -849,7 +848,11 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) colvarbias_restraint_moving::init(conf); colvarbias_restraint_k_moving::init(conf); - provide(f_cvb_scalar_variables); + get_keyval(conf, "lowerWallConstant", lower_wall_k, + (lower_wall_k > 0.0) ? lower_wall_k : force_k); + get_keyval(conf, "upperWallConstant", upper_wall_k, + (upper_wall_k > 0.0) ? upper_wall_k : force_k); + enable(f_cvb_scalar_variables); size_t i; @@ -857,9 +860,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) bool b_null_lower_walls = false; if (lower_walls.size() == 0) { b_null_lower_walls = true; - lower_walls.resize(number_of_colvars()); - for (i = 0; i < colvars.size(); i++) { - lower_walls[i].type(colvars[i]->value()); + lower_walls.resize(num_variables()); + for (i = 0; i < num_variables(); i++) { + lower_walls[i].type(variables(i)->value()); lower_walls[i].reset(); } } @@ -872,9 +875,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) bool b_null_upper_walls = false; if (upper_walls.size() == 0) { b_null_upper_walls = true; - upper_walls.resize(number_of_colvars()); - for (i = 0; i < colvars.size(); i++) { - upper_walls[i].type(colvars[i]->value()); + upper_walls.resize(num_variables()); + for (i = 0; i < num_variables(); i++) { + upper_walls[i].type(variables(i)->value()); upper_walls[i].reset(); } } @@ -890,17 +893,17 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) } if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) { - for (i = 0; i < colvars.size(); i++) { - if (colvars[i]->is_enabled(f_cv_periodic)) { + for (i = 0; i < num_variables(); i++) { + if (variables(i)->is_enabled(f_cv_periodic)) { cvm::error("Error: at least one variable is periodic, " - "both walls must be provided .\n", INPUT_ERROR); + "both walls must be provided.\n", INPUT_ERROR); return INPUT_ERROR; } } } if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) { - for (i = 0; i < colvars.size(); i++) { + for (i = 0; i < num_variables(); i++) { if (lower_walls[i] >= upper_walls[i]) { cvm::error("Error: one upper wall, "+ cvm::to_str(upper_walls[i])+ @@ -909,13 +912,24 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) INPUT_ERROR); } } + if (lower_wall_k * upper_wall_k == 0.0) { + cvm::error("Error: lowerWallConstant and upperWallConstant, " + "when defined, must both be positive.\n", + INPUT_ERROR); + return INPUT_ERROR; + } + force_k = lower_wall_k * upper_wall_k; + // transform the two constants to relative values + // (allow changing both via force_k) + lower_wall_k /= force_k; + upper_wall_k /= force_k; } - for (i = 0; i < colvars.size(); i++) { - if (colvars[i]->width != 1.0) - cvm::log("The force constant for colvar \""+colvars[i]->name+ + for (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 / (colvars[i]->width * colvars[i]->width))+ + cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+ " according to the specified width.\n"); } @@ -935,20 +949,20 @@ int colvarbias_restraint_harmonic_walls::update() void colvarbias_restraint_harmonic_walls::communicate_forces() { - for (size_t i = 0; i < colvars.size(); i++) { + for (size_t i = 0; i < num_variables(); i++) { if (cvm::debug()) { cvm::log("Communicating a force to colvar \""+ - colvars[i]->name+"\".\n"); + variables(i)->name+"\".\n"); } - colvars[i]->add_bias_force_actual_value(colvar_forces[i]); + variables(i)->add_bias_force_actual_value(colvar_forces[i]); } } cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const { - colvar *cv = colvars[i]; - colvarvalue const &cvv = colvars[i]->actual_value(); + colvar *cv = variables(i); + colvarvalue const &cvv = variables(i)->actual_value(); // For a periodic colvar, both walls may be applicable at the same time // in which case we pick the closer one @@ -958,21 +972,21 @@ cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const cvm::real const upper_wall_dist2 = cv->dist2(cvv, upper_walls[i]); if (lower_wall_dist2 < upper_wall_dist2) { cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); - if (grad < 0.0) { return grad; } + if (grad < 0.0) { return 0.5 * grad; } } else { cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); - if (grad > 0.0) { return grad; } + if (grad > 0.0) { return 0.5 * grad; } } return 0.0; } if (lower_walls.size() > 0) { cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); - if (grad < 0.0) { return grad; } + if (grad < 0.0) { return 0.5 * grad; } } if (upper_walls.size() > 0) { cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); - if (grad > 0.0) { return grad; } + if (grad > 0.0) { return 0.5 * grad; } } return 0.0; } @@ -981,7 +995,8 @@ cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) const { cvm::real const dist = colvar_distance(i); - return 0.5 * force_k / (colvars[i]->width * colvars[i]->width) * + cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k; + return 0.5 * force_k * scale / (variables(i)->width * variables(i)->width) * dist * dist; } @@ -989,15 +1004,16 @@ cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) con colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const { cvm::real const dist = colvar_distance(i); - return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) * - dist; + cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k; + return - force_k * scale / (variables(i)->width * variables(i)->width) * dist; } cvm::real colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i) const { cvm::real const dist = colvar_distance(i); - return 0.5 / (colvars[i]->width * colvars[i]->width) * + cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k; + return 0.5 * scale / (variables(i)->width * variables(i)->width) * dist * dist; } @@ -1054,16 +1070,16 @@ int colvarbias_restraint_linear::init(std::string const &conf) colvarbias_restraint_centers_moving::init(conf); colvarbias_restraint_k_moving::init(conf); - for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->is_enabled(f_cv_periodic)) { + for (size_t i = 0; i < num_variables(); i++) { + if (variables(i)->is_enabled(f_cv_periodic)) { cvm::error("Error: linear biases cannot be applied to periodic variables.\n", INPUT_ERROR); return INPUT_ERROR; } - if (colvars[i]->width != 1.0) - cvm::log("The force constant for colvar \""+colvars[i]->name+ + 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 / colvars[i]->width)+ + cvm::to_str(force_k / variables(i)->width)+ " according to the specified width.\n"); } @@ -1113,19 +1129,19 @@ cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const { - return force_k / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]); + return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]); } colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const { - return -1.0 * force_k / colvars[i]->width; + return -1.0 * force_k / variables(i)->width; } cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const { - return 1.0 / colvars[i]->width * (colvars[i]->value() - colvar_centers[i]); + return 1.0 / variables(i)->width * (variables(i)->value() - colvar_centers[i]); } @@ -1279,16 +1295,16 @@ int colvarbias_restraint_histogram::update() size_t vector_size = 0; size_t icv; - for (icv = 0; icv < colvars.size(); icv++) { - vector_size += colvars[icv]->value().size(); + for (icv = 0; icv < num_variables(); icv++) { + vector_size += variables(icv)->value().size(); } cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size); // calculate the histogram p.reset(); - for (icv = 0; icv < colvars.size(); icv++) { - colvarvalue const &cv = colvars[icv]->value(); + for (icv = 0; icv < num_variables(); icv++) { + colvarvalue const &cv = variables(icv)->value(); if (cv.type() == colvarvalue::type_scalar) { cvm::real const cv_value = cv.real_value; size_t igrid; @@ -1309,7 +1325,9 @@ int colvarbias_restraint_histogram::update() } } } else { - // TODO + cvm::error("Error: unsupported type for variable "+variables(icv)->name+".\n", + COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; } } @@ -1320,8 +1338,8 @@ int colvarbias_restraint_histogram::update() bias_energy = 0.5 * force_k_cv * p_diff * p_diff; // calculate the forces - for (icv = 0; icv < colvars.size(); icv++) { - colvarvalue const &cv = colvars[icv]->value(); + for (icv = 0; icv < num_variables(); icv++) { + colvarvalue const &cv = variables(icv)->value(); colvarvalue &cv_force = colvar_forces[icv]; cv_force.type(cv); cv_force.reset(); @@ -1363,10 +1381,10 @@ int colvarbias_restraint_histogram::update() std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) { if (b_write_histogram) { - std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat"); + std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat"); std::ofstream os(file_name.c_str()); - os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width) - << " " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3) + os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width) + << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3) << ")\n"; size_t igrid; for (igrid = 0; igrid < p.size(); igrid++) { diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index cbdd9c44d1..98b967abdb 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -260,6 +260,12 @@ protected: /// \brief Location of the upper walls std::vector upper_walls; + /// \brief If both walls are defined, use this k for the lower + cvm::real lower_wall_k; + + /// \brief If both walls are defined, use this k for the upper + cvm::real upper_wall_k; + virtual cvm::real colvar_distance(size_t i) const; virtual cvm::real restraint_potential(size_t i) const; virtual colvarvalue const restraint_force(size_t i) const; diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index d99da56386..786bc032d2 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -48,8 +48,6 @@ colvar::cvc::cvc(std::string const &conf) get_keyval(conf, "period", period, 0.0); get_keyval(conf, "wrapAround", wrap_center, 0.0); - // All cvcs implement this - provide(f_cvc_debug_gradient); get_keyval_feature((colvarparse *)this, conf, "debugGradients", f_cvc_debug_gradient, false, parse_silent); @@ -63,6 +61,8 @@ colvar::cvc::cvc(std::string const &conf) int colvar::cvc::init_total_force_params(std::string const &conf) { + if (cvm::get_error()) return COLVARS_ERROR; + if (get_keyval_feature(this, conf, "oneSiteSystemForce", f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) { cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: " @@ -72,6 +72,19 @@ int colvar::cvc::init_total_force_params(std::string const &conf) f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) { cvm::log("Computing total force on group 1 only"); } + + if (! is_enabled(f_cvc_one_site_total_force)) { + // check whether any of the other atom groups is dummy + std::vector::iterator agi = atom_groups.begin(); + agi++; + for ( ; agi != atom_groups.end(); agi++) { + if ((*agi)->b_dummy) { + provide(f_cvc_inv_gradient, false); + provide(f_cvc_Jacobian, false); + } + } + } + return COLVARS_OK; } @@ -87,8 +100,7 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, group->key = group_key; if (b_try_scalable) { - // TODO rewrite this logic in terms of dependencies - if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) { + if (is_available(f_cvc_scalable_com) && is_enabled(f_cvc_com_based)) { enable(f_cvc_scalable_com); enable(f_cvc_scalable); // The CVC makes the feature available; diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index a230cae8a9..ec215cbad1 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -343,6 +343,15 @@ public: virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force(colvarvalue const &force); + /// Redefined to override the distance ones + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to override the distance ones + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to override the distance ones + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; }; diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 4c3390a8bd..0204f3b4b1 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -21,7 +21,7 @@ colvar::angle::angle(std::string const &conf) function_type = "angle"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); @@ -40,7 +40,7 @@ colvar::angle::angle(cvm::atom const &a1, function_type = "angle"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); group1 = new cvm::atom_group(std::vector(1, a1)); group2 = new cvm::atom_group(std::vector(1, a2)); @@ -259,7 +259,7 @@ colvar::dihedral::dihedral(std::string const &conf) b_periodic = true; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); @@ -285,7 +285,7 @@ colvar::dihedral::dihedral(cvm::atom const &a1, b_periodic = true; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); b_1site_force = false; diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 7481dd360b..f46270246f 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -23,7 +23,7 @@ colvar::distance::distance(std::string const &conf) function_type = "distance"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); @@ -44,7 +44,7 @@ colvar::distance::distance() function_type = "distance"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); b_no_PBC = false; x.type(colvarvalue::type_scalar); } @@ -106,7 +106,7 @@ colvar::distance_vec::distance_vec(std::string const &conf) : distance(conf) { function_type = "distance_vec"; - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_3vector); } @@ -115,7 +115,7 @@ colvar::distance_vec::distance_vec() : distance() { function_type = "distance_vec"; - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_3vector); } @@ -176,7 +176,7 @@ colvar::distance_z::distance_z(std::string const &conf) function_type = "distance_z"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_scalar); // TODO detect PBC from MD engine (in simple cases) @@ -228,7 +228,7 @@ colvar::distance_z::distance_z() function_type = "distance_z"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_scalar); } @@ -372,7 +372,7 @@ colvar::distance_xy::distance_xy(std::string const &conf) function_type = "distance_xy"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_scalar); } @@ -383,7 +383,7 @@ colvar::distance_xy::distance_xy() function_type = "distance_xy"; provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_scalar); } @@ -479,7 +479,7 @@ colvar::distance_dir::distance_dir(std::string const &conf) : distance(conf) { function_type = "distance_dir"; - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_unit3vector); } @@ -488,7 +488,7 @@ colvar::distance_dir::distance_dir() : distance() { function_type = "distance_dir"; - provide(f_cvc_com_based); + enable(f_cvc_com_based); x.type(colvarvalue::type_unit3vector); } @@ -529,6 +529,27 @@ void colvar::distance_dir::apply_force(colvarvalue const &force) } +cvm::real colvar::distance_dir::dist2(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return (x1.rvector_value - x2.rvector_value).norm2(); +} + + +colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector); +} + + +colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector); +} + + colvar::distance_inv::distance_inv(std::string const &conf) : distance(conf) diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index e3ccd4ecf5..8252f77e62 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -14,9 +14,6 @@ colvardeps::~colvardeps() { size_t i; - for (i=0; iavailable = true; +void colvardeps::provide(int feature_id, bool truefalse) { + feature_states[feature_id].available = truefalse; +} + + +void colvardeps::set_enabled(int feature_id, bool truefalse) { +// if (!is_static(feature_id)) { +// cvm::error("Cannot set feature " + features()[feature_id]->description + " statically in " + description + ".\n"); +// return; +// } + if (truefalse) { + // Resolve dependencies too + enable(feature_id); + } else { + feature_states[feature_id].enabled = false; + } } bool colvardeps::get_keyval_feature(colvarparse *cvp, - std::string const &conf, char const *key, - int feature_id, bool const &def_value, - colvarparse::Parse_Mode const parse_mode) + std::string const &conf, char const *key, + int feature_id, bool const &def_value, + colvarparse::Parse_Mode const parse_mode) { + if (!is_user(feature_id)) { + cvm::error("Cannot set feature " + features()[feature_id]->description + " from user input in " + description + ".\n"); + return false; + } bool value; bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode); if (value) enable(feature_id); @@ -52,19 +67,19 @@ 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 + 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; bool ok; feature *f = features()[feature_id]; - feature_state *fs = feature_states[feature_id]; + feature_state *fs = &feature_states[feature_id]; if (cvm::debug()) { cvm::log("DEPS: " + description + @@ -88,6 +103,14 @@ int colvardeps::enable(int feature_id, return COLVARS_ERROR; } + if (!toplevel && !is_dynamic(feature_id)) { + if (!dry_run) { + cvm::log("Non-dynamic feature : \"" + f->description + + "\" in " + description + " may not be enabled as a dependency.\n"); + } + return COLVARS_ERROR; + } + // 1) enforce exclusions for (i=0; irequires_exclude.size(); i++) { feature *g = features()[f->requires_exclude[i]]; @@ -168,9 +191,9 @@ int colvardeps::enable(int feature_id, if (res != COLVARS_OK) { if (!dry_run) { cvm::log("...required by \"" + f->description + "\" in " + description); - } - if (toplevel) { - cvm::error("Error: Failed dependency in " + description + "."); + if (toplevel) { + cvm::error("Error: Failed dependency in " + description + "."); + } } return res; } @@ -194,9 +217,12 @@ int colvardeps::enable(int feature_id, // // we need refs to parents to walk up the deps tree! // // or refresh // } +void colvardeps::init_feature(int feature_id, const char *description, feature_type type) { + features()[feature_id]->description = description; + features()[feature_id]->type = type; +} - // Shorthand macros for describing dependencies -#define f_description(f, d) features()[f]->description = d +// 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); \ @@ -216,35 +242,31 @@ void colvardeps::init_cvb_requires() { for (i = 0; i < f_cvb_ntot; i++) { 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_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_req_children(f_cvb_get_total_force, f_cv_total_force); + + init_feature(f_cvb_history_dependent, "history-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); } - f_description(f_cvb_active, "active"); - f_req_children(f_cvb_active, f_cv_active); - - f_description(f_cvb_apply_force, "apply force"); - f_req_children(f_cvb_apply_force, f_cv_gradient); - - f_description(f_cvb_get_total_force, "obtain total force"); - f_req_children(f_cvb_get_total_force, f_cv_total_force); - - f_description(f_cvb_history_dependent, "history-dependent"); - - f_description(f_cvb_scalar_variables, "require scalar variables"); - f_req_children(f_cvb_scalar_variables, f_cv_scalar); - // Initialize feature_states for each instance feature_states.reserve(f_cvb_ntot); for (i = 0; i < f_cvb_ntot; i++) { - feature_states.push_back(new feature_state(true, false)); + feature_states.push_back(feature_state(true, false)); // Most features are available, so we set them so // and list exceptions below } - - // some biases are not history-dependent - feature_states[f_cvb_history_dependent]->available = false; - - // by default, biases should work with vector variables, too - feature_states[f_cvb_scalar_variables]->available = false; } @@ -255,117 +277,111 @@ void colvardeps::init_cv_requires() { features().push_back(new feature); } - f_description(f_cv_active, "active"); + init_feature(f_cv_active, "active", f_type_dynamic); f_req_children(f_cv_active, f_cvc_active); // Colvars must be either a linear combination, or scalar (and polynomial) or scripted f_req_alt3(f_cv_active, f_cv_scalar, f_cv_linear, f_cv_scripted); - f_description(f_cv_gradient, "gradient"); + init_feature(f_cv_gradient, "gradient", f_type_dynamic); f_req_children(f_cv_gradient, f_cvc_gradient); - f_description(f_cv_collect_gradient, "collect 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); - f_description(f_cv_fdiff_velocity, "fdiff_velocity"); + init_feature(f_cv_fdiff_velocity, "fdiff_velocity", f_type_dynamic); // System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly - f_description(f_cv_total_force, "total force"); + 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 - f_description(f_cv_total_force_calc, "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); - f_description(f_cv_Jacobian, "Jacobian derivative"); + 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); - f_description(f_cv_hide_Jacobian, "hide Jacobian force"); + 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 - f_description(f_cv_extended_Lagrangian, "extended Lagrangian"); + 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); - f_description(f_cv_Langevin, "Langevin dynamics"); + init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user); f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian); - f_description(f_cv_linear, "linear"); + init_feature(f_cv_linear, "linear", f_type_static); - f_description(f_cv_scalar, "scalar"); + init_feature(f_cv_scalar, "scalar", f_type_static); - f_description(f_cv_output_energy, "output energy"); + init_feature(f_cv_output_energy, "output energy", f_type_user); - f_description(f_cv_output_value, "output value"); + init_feature(f_cv_output_value, "output value", f_type_user); - f_description(f_cv_output_velocity, "output velocity"); + init_feature(f_cv_output_velocity, "output velocity", f_type_user); f_req_self(f_cv_output_velocity, f_cv_fdiff_velocity); - f_description(f_cv_output_applied_force, "output applied force"); + init_feature(f_cv_output_applied_force, "output applied force", f_type_user); - f_description(f_cv_output_total_force, "output total force"); + 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); - f_description(f_cv_subtract_applied_force, "subtract applied force from 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); - f_description(f_cv_lower_boundary, "lower boundary"); + init_feature(f_cv_lower_boundary, "lower boundary", f_type_user); f_req_self(f_cv_lower_boundary, f_cv_scalar); - f_description(f_cv_upper_boundary, "upper boundary"); + init_feature(f_cv_upper_boundary, "upper boundary", f_type_user); f_req_self(f_cv_upper_boundary, f_cv_scalar); - f_description(f_cv_grid, "grid"); + init_feature(f_cv_grid, "grid", f_type_user); f_req_self(f_cv_grid, f_cv_lower_boundary); f_req_self(f_cv_grid, f_cv_upper_boundary); - f_description(f_cv_lower_wall, "lower wall"); - f_req_self(f_cv_lower_wall, f_cv_lower_boundary); - f_req_self(f_cv_lower_wall, f_cv_gradient); + init_feature(f_cv_runave, "running average", f_type_user); - f_description(f_cv_upper_wall, "upper wall"); - f_req_self(f_cv_upper_wall, f_cv_upper_boundary); - f_req_self(f_cv_upper_wall, f_cv_gradient); + init_feature(f_cv_corrfunc, "correlation function", f_type_user); - f_description(f_cv_runave, "running average"); - - f_description(f_cv_corrfunc, "correlation function"); - - // The features below are set programmatically - f_description(f_cv_scripted, "scripted"); - f_description(f_cv_periodic, "periodic"); + init_feature(f_cv_scripted, "scripted", f_type_static); + init_feature(f_cv_periodic, "periodic", f_type_static); f_req_self(f_cv_periodic, f_cv_homogeneous); - f_description(f_cv_scalar, "scalar"); - f_description(f_cv_linear, "linear"); - f_description(f_cv_homogeneous, "homogeneous"); + 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); } // Initialize feature_states for each instance feature_states.reserve(f_cv_ntot); for (i = 0; i < f_cv_ntot; i++) { - feature_states.push_back(new feature_state(true, false)); + feature_states.push_back(feature_state(true, false)); // Most features are available, so we set them so // and list exceptions below } - // properties that may NOT be enabled as a dependency - int unavailable_deps[] = { - f_cv_lower_boundary, - f_cv_upper_boundary, - f_cv_extended_Lagrangian, - f_cv_Langevin, - f_cv_scripted, - f_cv_periodic, - f_cv_scalar, - f_cv_linear, - f_cv_homogeneous - }; - for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) { - feature_states[unavailable_deps[i]]->available = false; - } +// // properties that may NOT be enabled as a dependency +// // This will be deprecated by feature types +// int unavailable_deps[] = { +// f_cv_lower_boundary, +// f_cv_upper_boundary, +// f_cv_extended_Lagrangian, +// f_cv_Langevin, +// f_cv_scripted, +// f_cv_periodic, +// f_cv_scalar, +// f_cv_linear, +// f_cv_homogeneous +// }; +// for (i = 0; i < sizeof(unavailable_deps) / sizeof(unavailable_deps[0]); i++) { +// feature_states[unavailable_deps[i]].available = false; +// } } @@ -377,34 +393,34 @@ void colvardeps::init_cvc_requires() { features().push_back(new feature); } - f_description(f_cvc_active, "active"); + 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); - f_description(f_cvc_scalar, "scalar"); + init_feature(f_cvc_scalar, "scalar", f_type_static); - f_description(f_cvc_gradient, "gradient"); + init_feature(f_cvc_gradient, "gradient", f_type_dynamic); - f_description(f_cvc_inv_gradient, "inverse gradient"); + init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic); f_req_self(f_cvc_inv_gradient, f_cvc_gradient); - f_description(f_cvc_debug_gradient, "debug gradient"); + init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user); f_req_self(f_cvc_debug_gradient, f_cvc_gradient); - f_description(f_cvc_Jacobian, "Jacobian derivative"); + init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic); f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient); - f_description(f_cvc_com_based, "depends on group centers of mass"); + init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static); // Compute total force on first site only to avoid unwanted // coupling to other colvars (see e.g. Ciccotti et al., 2005) - f_description(f_cvc_one_site_total_force, "compute total collective force only from one group center"); + init_feature(f_cvc_one_site_total_force, "compute total collective force only from one group center", f_type_user); f_req_self(f_cvc_one_site_total_force, f_cvc_com_based); - f_description(f_cvc_scalable, "scalable calculation"); + init_feature(f_cvc_scalable, "scalable calculation", f_type_static); f_req_self(f_cvc_scalable, f_cvc_scalable_com); - f_description(f_cvc_scalable_com, "scalable calculation of centers of mass"); + 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); @@ -414,23 +430,25 @@ void colvardeps::init_cvc_requires() { } // Initialize feature_states for each instance - // default as unavailable, not enabled + // 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++) { - feature_states.push_back(new feature_state(false, false)); + 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; + feature_states[f_cvc_active].available = true; + feature_states[f_cvc_gradient].available = true; // Features that are implemented by default if their requirements are - feature_states[f_cvc_one_site_total_force]->available = true; + 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; + 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; } @@ -442,21 +460,21 @@ void colvardeps::init_ag_requires() { features().push_back(new feature); } - f_description(f_ag_active, "active"); - f_description(f_ag_center, "translational fit"); - f_description(f_ag_rotate, "rotational fit"); - f_description(f_ag_fitting_group, "reference positions group"); - f_description(f_ag_fit_gradient_group, "fit gradient for main group"); - f_description(f_ag_fit_gradient_ref, "fit gradient for reference group"); - f_description(f_ag_atom_forces, "atomic forces"); + 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_fit_gradient_group, "fit gradient for main group", f_type_static); + init_feature(f_ag_fit_gradient_ref, "fit gradient for reference group", f_type_static); + 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 - f_description(f_ag_scalable, "scalable group calculation"); - f_description(f_ag_scalable_com, "scalable group center of mass calculation"); + 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); -// f_description(f_ag_min_msd_fit, "minimum MSD fit") +// 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) @@ -466,15 +484,15 @@ void colvardeps::init_ag_requires() { // default as unavailable, not enabled feature_states.reserve(f_ag_ntot); for (i = 0; i < colvardeps::f_ag_ntot; i++) { - feature_states.push_back(new feature_state(false, false)); + 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; + 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; + 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_scalable].available = true; } @@ -482,7 +500,7 @@ void colvardeps::print_state() { size_t i; cvm::log("Enabled features of " + description); for (i = 0; i < feature_states.size(); i++) { - if (feature_states[i]->enabled) + if (feature_states[i].enabled) cvm::log("- " + features()[i]->description); } for (i=0; i feature_states; + +private: + /// List of the states of all features + std::vector feature_states; + + /// Enum of possible feature types + enum feature_type { + f_type_not_set, + f_type_dynamic, + f_type_user, + f_type_static + }; + +public: + /// 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); /// Describes a feature and its dependecies /// used in a static array within each subclass @@ -83,8 +91,18 @@ public: // features that this feature requires in children std::vector requires_children; + + inline bool is_dynamic() { return type == f_type_dynamic; } + inline bool is_static() { return type == f_type_static; } + inline bool is_user() { return type == f_type_user; } + /// Type of this feature, from the enum feature_type + feature_type type; }; + 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; } + // Accessor to array of all features with deps, static in most derived classes // Subclasses with dynamic dependency trees may override this // with a non-static array @@ -100,9 +118,8 @@ public: /// (useful for cvcs because their children are member objects) void remove_all_children(); - - private: + // pointers to objects this object depends on // list should be maintained by any code that modifies the object // this could be secured by making lists of colvars / cvcs / atom groups private and modified through accessor functions @@ -130,16 +147,28 @@ public: // Checks whether given feature is enabled // Defaults to querying f_*_active inline bool is_enabled(int f = f_cv_active) const { - return feature_states[f]->enabled; + return feature_states[f].enabled; } // Checks whether given feature is available // Defaults to querying f_*_active inline bool is_available(int f = f_cv_active) const { - return feature_states[f]->available; + return feature_states[f].available; } - void provide(int feature_id); // set the feature's flag to available in local object + /// Set the feature's available flag, without checking + /// To be used for dynamic properties + /// dependencies will be checked by enable() + void provide(int feature_id, bool truefalse = true); + + /// Set the feature's enabled flag, without dependency check or resolution + /// To be used for static properties only + /// Checking for availability is up to the caller + void set_enabled(int feature_id, bool truefalse = true); + +protected: + + /// Parse a keyword and enable a feature accordingly bool get_keyval_feature(colvarparse *cvp, @@ -147,6 +176,8 @@ public: int feature_id, bool const &def_value, colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal); +public: + int enable(int f, bool dry_run = false, bool toplevel = true); // enable a feature and recursively solve its dependencies // dry_run is set to true to recursively test if a feature is available, without enabling it // int disable(int f); @@ -165,6 +196,7 @@ public: f_cvb_get_total_force, // requires total forces f_cvb_history_dependent, // depends on simulation history f_cvb_scalar_variables, // requires scalar colvars + f_cvb_calc_pmf, // whether this bias will compute a PMF f_cvb_ntot }; @@ -216,12 +248,6 @@ public: /// be used by the biases or in analysis (needs lower and upper /// boundary) f_cv_grid, - /// \brief Apply a restraining potential (|x-xb|^2) to the colvar - /// when it goes below the lower wall - f_cv_lower_wall, - /// \brief Apply a restraining potential (|x-xb|^2) to the colvar - /// when it goes above the upper wall - f_cv_upper_wall, /// \brief Compute running average f_cv_runave, /// \brief Compute time correlation function diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index b9a435152b..10cd3c0e47 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -36,6 +36,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) "variable module twice.\n"); return; } + + depth_s = 0; + cvm::log(cvm::line_marker); cvm::log("Initializing the collective variables module, version "+ cvm::to_str(COLVARS_VERSION)+".\n"); @@ -71,6 +74,48 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) } +colvarmodule * colvarmodule::main() +{ + return proxy->colvars; +} + + +std::vector *colvarmodule::variables() +{ + return &colvars; +} + + +std::vector *colvarmodule::variables_active() +{ + return &colvars_active; +} + + +std::vector *colvarmodule::variables_active_smp() +{ + return &colvars_smp; +} + + +std::vector *colvarmodule::variables_active_smp_items() +{ + return &colvars_smp_items; +} + + +std::vector *colvarmodule::biases_active() +{ + return &(biases_active_); +} + + +size_t colvarmodule::size() const +{ + return colvars.size() + biases.size(); +} + + int colvarmodule::read_config_file(char const *config_filename) { cvm::log(cvm::line_marker); @@ -118,8 +163,29 @@ int colvarmodule::read_config_string(std::string const &config_str) } +std::istream & colvarmodule::getline(std::istream &is, std::string &line) +{ + std::string l; + if (std::getline(is, l)) { + size_t const sz = l.size(); + if (sz > 0) { + if (l[sz-1] == '\r' ) { + line = l.substr(0, sz-1); + } else { + line = l; + } + } else { + line.clear(); + } + } + return is; +} + + int colvarmodule::parse_config(std::string &conf) { + extra_conf.clear(); + // parse global options if (catch_input_errors(parse_global_params(conf))) { return get_error(); @@ -140,6 +206,15 @@ int colvarmodule::parse_config(std::string &conf) return get_error(); } + if (extra_conf.size()) { + catch_input_errors(parse_global_params(extra_conf)); + catch_input_errors(parse_colvars(extra_conf)); + catch_input_errors(parse_biases(extra_conf)); + parse->check_keywords(extra_conf, "colvarmodule"); + extra_conf.clear(); + if (get_error() != COLVARS_OK) return get_error(); + } + cvm::log(cvm::line_marker); cvm::log("Collective variables module (re)initialized.\n"); cvm::log(cvm::line_marker); @@ -156,11 +231,20 @@ int colvarmodule::parse_config(std::string &conf) } +int colvarmodule::append_new_config(std::string const &new_conf) +{ + extra_conf += new_conf; + return COLVARS_OK; +} + + int colvarmodule::parse_global_params(std::string const &conf) { + colvarmodule *cvm = cvm::main(); + std::string index_file_name; if (parse->get_keyval(conf, "indexFile", index_file_name)) { - read_index_file(index_file_name.c_str()); + cvm->read_index_file(index_file_name.c_str()); } if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) { @@ -216,8 +300,8 @@ int colvarmodule::parse_colvars(std::string const &conf) if (colvar_conf.size()) { cvm::log(cvm::line_marker); cvm::increase_depth(); - colvars.push_back(new colvar(colvar_conf)); - if (cvm::get_error() || + colvars.push_back(new colvar()); + if (((colvars.back())->init(colvar_conf) != COLVARS_OK) || ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) { cvm::log("Error while constructing colvar number " + cvm::to_str(colvars.size()) + " : deleting."); @@ -262,8 +346,7 @@ bool colvarmodule::check_new_bias(std::string &conf, char const *key) template int colvarmodule::parse_biases_type(std::string const &conf, - char const *keyword, - size_t &bias_count) + char const *keyword) { std::string bias_conf = ""; size_t conf_saved_pos = 0; @@ -277,7 +360,6 @@ int colvarmodule::parse_biases_type(std::string const &conf, return COLVARS_ERROR; } cvm::decrease_depth(); - bias_count++; } else { cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n", INPUT_ERROR); @@ -295,28 +377,28 @@ int colvarmodule::parse_biases(std::string const &conf) cvm::log("Initializing the collective variables biases.\n"); /// initialize ABF instances - parse_biases_type(conf, "abf", n_abf_biases); + parse_biases_type(conf, "abf"); /// initialize adaptive linear biases - parse_biases_type(conf, "ALB", n_rest_biases); + parse_biases_type(conf, "ALB"); /// initialize harmonic restraints - parse_biases_type(conf, "harmonic", n_rest_biases); + parse_biases_type(conf, "harmonic"); /// initialize harmonic walls restraints - parse_biases_type(conf, "harmonicWalls", n_rest_biases); + parse_biases_type(conf, "harmonicWalls"); /// initialize histograms - parse_biases_type(conf, "histogram", n_histo_biases); + parse_biases_type(conf, "histogram"); /// initialize histogram restraints - parse_biases_type(conf, "histogramRestraint", n_rest_biases); + parse_biases_type(conf, "histogramRestraint"); /// initialize linear restraints - parse_biases_type(conf, "linear", n_rest_biases); + parse_biases_type(conf, "linear"); /// initialize metadynamics instances - parse_biases_type(conf, "metadynamics", n_meta_biases); + parse_biases_type(conf, "metadynamics"); if (use_scripted_forces) { cvm::log(cvm::line_marker); @@ -363,6 +445,36 @@ int colvarmodule::parse_biases(std::string const &conf) } +int colvarmodule::num_biases_feature(int feature_id) +{ + colvarmodule *cv = cvm::main(); + size_t n = 0; + for (std::vector::iterator bi = cv->biases.begin(); + bi != cv->biases.end(); + bi++) { + if ((*bi)->is_enabled(feature_id)) { + n++; + } + } + return n; +} + + +int colvarmodule::num_biases_type(std::string const &type) +{ + colvarmodule *cv = cvm::main(); + size_t n = 0; + for (std::vector::iterator bi = cv->biases.begin(); + bi != cv->biases.end(); + bi++) { + if ((*bi)->bias_type == type) { + n++; + } + } + return n; +} + + int colvarmodule::catch_input_errors(int result) { if (result != COLVARS_OK || get_error()) { @@ -376,8 +488,9 @@ int colvarmodule::catch_input_errors(int result) colvarbias * colvarmodule::bias_by_name(std::string const &name) { - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); + colvarmodule *cv = cvm::main(); + for (std::vector::iterator bi = cv->biases.begin(); + bi != cv->biases.end(); bi++) { if ((*bi)->name == name) { return (*bi); @@ -388,8 +501,9 @@ colvarbias * colvarmodule::bias_by_name(std::string const &name) { colvar *colvarmodule::colvar_by_name(std::string const &name) { - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); + colvarmodule *cv = cvm::main(); + for (std::vector::iterator cvi = cv->colvars.begin(); + cvi != cv->colvars.end(); cvi++) { if ((*cvi)->name == name) { return (*cvi); @@ -555,9 +669,15 @@ int colvarmodule::calc_colvars() int error_code = COLVARS_OK; std::vector::iterator cvi; - // Determine which colvars are active at this time step - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->feature_states[colvardeps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 0); + // Determine which colvars are active at this iteration + variables_active()->resize(0); + variables_active()->reserve(variables()->size()); + for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) { + // This is a dynamic feature - the next call should be to enable() + // or disable() when dynamic dependency resolution is fully implemented + (*cvi)->set_enabled(colvardeps::f_cv_active, + step_absolute() % (*cvi)->get_time_step_factor() == 0); + variables_active()->push_back(*cvi); } // if SMP support is available, split up the work @@ -565,25 +685,24 @@ int colvarmodule::calc_colvars() // first, calculate how much work (currently, how many active CVCs) each colvar has - colvars_smp.resize(0); - colvars_smp_items.resize(0); + variables_active_smp()->resize(0); + variables_active_smp_items()->resize(0); - colvars_smp.reserve(colvars.size()); - colvars_smp_items.reserve(colvars.size()); + variables_active_smp()->reserve(variables_active()->size()); + variables_active_smp_items()->reserve(variables_active()->size()); // set up a vector containing all components cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) { - if (!(*cvi)->is_enabled()) continue; error_code |= (*cvi)->update_cvc_flags(); size_t num_items = (*cvi)->num_active_cvcs(); - colvars_smp.reserve(colvars_smp.size() + num_items); - colvars_smp_items.reserve(colvars_smp_items.size() + num_items); + variables_active_smp()->reserve(variables_active_smp()->size() + num_items); + variables_active_smp_items()->reserve(variables_active_smp_items()->size() + num_items); for (size_t icvc = 0; icvc < num_items; icvc++) { - colvars_smp.push_back(*cvi); - colvars_smp_items.push_back(icvc); + variables_active_smp()->push_back(*cvi); + variables_active_smp_items()->push_back(icvc); } } cvm::decrease_depth(); @@ -592,8 +711,7 @@ int colvarmodule::calc_colvars() error_code |= proxy->smp_colvars_loop(); cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - if (!(*cvi)->is_enabled()) continue; + for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) { error_code |= (*cvi)->collect_cvc_data(); } cvm::decrease_depth(); @@ -602,8 +720,7 @@ int colvarmodule::calc_colvars() // calculate colvars one at a time cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - if (!(*cvi)->is_enabled()) continue; + for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) { error_code |= (*cvi)->calc(); if (cvm::get_error()) { return COLVARS_ERROR; @@ -627,6 +744,18 @@ int colvarmodule::calc_biases() std::vector::iterator bi; int error_code = COLVARS_OK; + // Total bias energy is reset before calling scripted biases + total_bias_energy = 0.0; + + // update the list of active biases + biases_active()->resize(0); + biases_active()->reserve(biases.size()); + for (bi = biases.begin(); bi != biases.end(); bi++) { + if ((*bi)->is_enabled()) { + biases_active()->push_back(*bi); + } + } + // if SMP support is available, split up the work if (proxy->smp_enabled() == COLVARS_OK) { @@ -645,7 +774,7 @@ int colvarmodule::calc_biases() } cvm::increase_depth(); - for (bi = biases.begin(); bi != biases.end(); bi++) { + for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) { error_code |= (*bi)->update(); if (cvm::get_error()) { return error_code; @@ -654,12 +783,10 @@ int colvarmodule::calc_biases() cvm::decrease_depth(); } - cvm::real total_bias_energy = 0.0; - for (bi = biases.begin(); bi != biases.end(); bi++) { + for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) { total_bias_energy += (*bi)->get_energy(); } - proxy->add_energy(total_bias_energy); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -675,7 +802,7 @@ int colvarmodule::update_colvar_forces() if (cvm::debug() && biases.size()) cvm::log("Collecting forces from all biases.\n"); cvm::increase_depth(); - for (bi = biases.begin(); bi != biases.end(); bi++) { + for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) { (*bi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; @@ -687,6 +814,11 @@ int colvarmodule::update_colvar_forces() error_code |= calc_scripted_forces(); } + // Now we have collected energies from both built-in and scripted biases + if (cvm::debug()) + cvm::log("Adding total bias energy: " + cvm::to_str(total_bias_energy) + "\n"); + proxy->add_energy(total_bias_energy); + cvm::real total_colvar_energy = 0.0; // sum up the forces for each colvar, including wall forces // and integrate any internal @@ -695,7 +827,7 @@ int colvarmodule::update_colvar_forces() cvm::log("Updating the internal degrees of freedom " "of colvars (if they have any).\n"); cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) { // Here we call even inactive colvars, so they accumulate biasing forces // as well as update their extended-system dynamics total_colvar_energy += (*cvi)->update_forces_energy(); @@ -704,6 +836,8 @@ int colvarmodule::update_colvar_forces() } } cvm::decrease_depth(); + if (cvm::debug()) + cvm::log("Adding total colvar energy: " + cvm::to_str(total_colvar_energy) + "\n"); proxy->add_energy(total_colvar_energy); // make collective variables communicate their forces to their @@ -711,9 +845,8 @@ int colvarmodule::update_colvar_forces() if (cvm::debug()) cvm::log("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); - for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) { if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) { - if (!(*cvi)->is_enabled()) continue; (*cvi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; @@ -800,8 +933,8 @@ int colvarmodule::analyze() cvm::log("Performing analysis.\n"); // perform colvar-specific analysis - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); + for (std::vector::iterator cvi = variables_active()->begin(); + cvi != variables_active()->end(); cvi++) { cvm::increase_depth(); (*cvi)->analyze(); @@ -825,8 +958,8 @@ 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 = colvars.begin(); - cvi != colvars.end(); cvi++) { + for (std::vector::iterator cvi = variables()->begin(); + cvi != variables()->end(); cvi++) { (*cvi)->setup(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -934,18 +1067,18 @@ int colvarmodule::setup_output() cvm::log("The restart output state file will be \""+restart_out_name+"\".\n"); } - output_prefix = proxy->output_prefix(); - if (output_prefix.size()) { + output_prefix() = proxy->output_prefix(); + if (output_prefix().size()) { cvm::log("The final output state file will be \""+ - (output_prefix.size() ? - std::string(output_prefix+".colvars.state") : + (output_prefix().size() ? + std::string(output_prefix()+".colvars.state") : std::string("colvars.state"))+"\".\n"); // cvm::log (cvm::line_marker); } cv_traj_name = - (output_prefix.size() ? - std::string(output_prefix+".colvars.traj") : + (output_prefix().size() ? + std::string(output_prefix()+".colvars.traj") : std::string("")); if (cv_traj_freq && cv_traj_name.size()) { @@ -1077,14 +1210,14 @@ continue the previous simulation.\n\n"); cvm::log(cvm::line_marker); // update this ahead of time in this special case - output_prefix = proxy->input_prefix(); - cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n"); + output_prefix() = proxy->input_prefix(); + cvm::log("All output files will now be saved with the prefix \""+output_prefix()+".tmp.*\".\n"); cvm::log(cvm::line_marker); cvm::log("Please review the important warning above. After that, you may rename:\n\ -\""+output_prefix+".tmp.colvars.state\"\n\ +\""+output_prefix()+".tmp.colvars.state\"\n\ to:\n\ \""+ proxy->input_prefix()+".colvars.state\"\n"); - output_prefix = output_prefix+".tmp"; + output_prefix() = output_prefix()+".tmp"; write_output_files(); cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR); } @@ -1104,8 +1237,8 @@ int colvarmodule::write_output_files() // if this is a simulation run (i.e. not a postprocessing), output data // must be written to be able to restart the simulation std::string const out_name = - (output_prefix.size() ? - std::string(output_prefix+".colvars.state") : + (output_prefix().size() ? + std::string(output_prefix()+".colvars.state") : std::string("colvars.state")); cvm::log("Saving collective variables state to \""+out_name+"\".\n"); @@ -1340,7 +1473,8 @@ std::ostream & colvarmodule::write_traj(std::ostream &os) void cvm::log(std::string const &message) { - size_t const d = depth(); + // allow logging when the module is not fully initialized + size_t const d = (cvm::main() != NULL) ? depth() : 0; if (d > 0) proxy->log((std::string(2*d, ' '))+message); else @@ -1365,19 +1499,20 @@ void cvm::decrease_depth() size_t & cvm::depth() { // NOTE: do not call log() or error() here, to avoid recursion - size_t const nt = proxy->smp_num_threads(); + colvarmodule *cv = cvm::main(); if (proxy->smp_enabled() == COLVARS_OK) { - if (depth_v.size() != nt) { - // update array of depths + int const nt = proxy->smp_num_threads(); + if (int(cv->depth_v.size()) != nt) { proxy->smp_lock(); - if (depth_v.size() > 0) { depth_s = depth_v[0]; } - depth_v.clear(); - depth_v.assign(nt, depth_s); + // update array of depths + if (cv->depth_v.size() > 0) { cv->depth_s = cv->depth_v[0]; } + cv->depth_v.clear(); + cv->depth_v.assign(nt, cv->depth_s); proxy->smp_unlock(); } - return depth_v[proxy->smp_thread_id()]; + return cv->depth_v[proxy->smp_thread_id()]; } - return depth_s; + return cv->depth_s; } @@ -1451,8 +1586,8 @@ int cvm::read_index_file(char const *filename) FILE_ERROR); } } - cvm::index_group_names.push_back(group_name); - cvm::index_groups.push_back(std::vector ()); + index_group_names.push_back(group_name); + index_groups.push_back(std::vector()); } else { cvm::error("Error: in parsing index file \""+ std::string(filename)+"\".\n", @@ -1462,7 +1597,7 @@ int cvm::read_index_file(char const *filename) int atom_number = 1; size_t pos = is.tellg(); while ( (is >> atom_number) && (atom_number > 0) ) { - (cvm::index_groups.back()).push_back(atom_number); + (index_groups.back()).push_back(atom_number); pos = is.tellg(); } is.clear(); @@ -1532,8 +1667,8 @@ int cvm::load_coords_xyz(char const *filename, + std::string(filename) + ".\n", INPUT_ERROR); } // skip comment line - std::getline(xyz_is, line); - std::getline(xyz_is, line); + cvm::getline(xyz_is, line); + cvm::getline(xyz_is, line); xyz_is.width(255); std::vector::iterator pos_i = pos.begin(); @@ -1543,7 +1678,7 @@ int cvm::load_coords_xyz(char const *filename, for ( ; pos_i != pos.end() ; pos_i++, index++) { while ( next < *index ) { - std::getline(xyz_is, line); + cvm::getline(xyz_is, line); next++; } xyz_is >> symbol; @@ -1558,15 +1693,9 @@ int cvm::load_coords_xyz(char const *filename, return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -// static pointers -std::vector colvarmodule::colvars; -std::vector colvarmodule::biases; -size_t colvarmodule::n_abf_biases = 0; -size_t colvarmodule::n_rest_biases = 0; -size_t colvarmodule::n_histo_biases = 0; -size_t colvarmodule::n_meta_biases = 0; -colvarproxy *colvarmodule::proxy = NULL; +// shared pointer to the proxy object +colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07; @@ -1575,14 +1704,9 @@ long colvarmodule::it = 0; long colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; -size_t colvarmodule::depth_s = 0; -std::vector colvarmodule::depth_v(0); bool colvarmodule::b_analysis = false; -std::list colvarmodule::index_group_names; -std::list > colvarmodule::index_groups; -bool colvarmodule::use_scripted_forces = false; -bool colvarmodule::scripting_after_biases = true; -std::string colvarmodule::output_prefix = ""; +bool colvarmodule::use_scripted_forces = false; +bool colvarmodule::scripting_after_biases = true; // i/o constants size_t const colvarmodule::it_width = 12; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 0e98709adb..b4f80e0b5e 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -11,7 +11,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2016-12-27" +#define COLVARS_VERSION "2017-03-09" #endif #ifndef COLVARS_DEBUG @@ -161,12 +161,37 @@ public: /// dt) static real debug_gradients_step_size; - /// Prefix for all output files for this run - static std::string output_prefix; +private: + /// Prefix for all output files for this run + std::string cvm_output_prefix; + +public: + /// Accessor for the above + static inline std::string &output_prefix() + { + colvarmodule *cv = colvarmodule::main(); + return cv->cvm_output_prefix; + } + +private: /// Array of collective variables - static std::vector colvars; + std::vector colvars; + + /// Array of collective variables + std::vector colvars_active; + + /// Collective variables to be calculated on different threads; + /// colvars with multple items (e.g. multiple active CVCs) are duplicated + std::vector colvars_smp; + /// Indexes of the items to calculate for each colvar + std::vector colvars_smp_items; + +public: + + /// Array of collective variables + std::vector *variables(); /* TODO: implement named CVCs /// Array of named (reusable) collective variable components @@ -177,26 +202,31 @@ public: } */ + /// Collective variables with the active flag on + std::vector *variables_active(); + /// Collective variables to be calculated on different threads; /// colvars with multple items (e.g. multiple active CVCs) are duplicated - std::vector colvars_smp; + std::vector *variables_active_smp(); + /// Indexes of the items to calculate for each colvar - std::vector colvars_smp_items; + std::vector *variables_active_smp_items(); /// Array of collective variable biases - static std::vector biases; - /// \brief Number of ABF biases initialized (in normal conditions - /// should be 1) - static size_t n_abf_biases; - /// \brief Number of metadynamics biases initialized (in normal - /// conditions should be 1) - static size_t n_meta_biases; - /// \brief Number of restraint biases initialized (no limit on the - /// number) - static size_t n_rest_biases; - /// \brief Number of histograms initialized (no limit on the - /// number) - static size_t n_histo_biases; + std::vector biases; + + /// Energy of built-in and scripted biases, summed per time-step + real total_bias_energy; + +private: + + /// Array of active collective variable biases + std::vector biases_active_; + +public: + + /// Array of active collective variable biases + std::vector *biases_active(); /// \brief Whether debug output should be enabled (compile-time option) static inline bool debug() @@ -205,10 +235,7 @@ public: } /// \brief How many objects are configured yet? - inline size_t size() const - { - return colvars.size() + biases.size(); - } + size_t size() const; /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name @@ -230,9 +257,12 @@ public: /// \brief Parse a "clean" config string (no comments) int parse_config(std::string &conf); - // Parse functions (setup internal data based on a string) + /// Allow reading from Windows text files using using std::getline + /// (which can still be used when the text is produced by Colvars itself) + static std::istream & getline(std::istream &is, std::string &line); + /// Parse the few module's global parameters int parse_global_params(std::string const &conf); @@ -242,14 +272,33 @@ public: /// Parse and initialize collective variable biases int parse_biases(std::string const &conf); + /// \brief Add new configuration during parsing (e.g. to implement + /// back-compatibility); cannot be nested, i.e. conf should not contain + /// anything that triggers another call + int append_new_config(std::string const &conf); + +private: + + /// Auto-generated configuration during parsing (e.g. to implement + /// back-compatibility) + std::string extra_conf; + /// Parse and initialize collective variable biases of a specific type template - int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count); + int parse_biases_type(std::string const &conf, char const *keyword); /// Test error condition and keyword parsing /// on error, delete new bias bool check_new_bias(std::string &conf, char const *key); +public: + + /// Return how many biases have this feature enabled + static int num_biases_feature(int feature_id); + + /// Return how many biases are defined with this type + static int num_biases_type(std::string const &type); + private: /// Useful wrapper to interrupt parsing if any error occurs int catch_input_errors(int result); @@ -449,13 +498,13 @@ public: /// \brief Names of groups from a Gromacs .ndx file to be read at startup - static std::list index_group_names; + std::list index_group_names; /// \brief Groups from a Gromacs .ndx file read at startup - static std::list > index_groups; + std::list > index_groups; /// \brief Read a Gromacs .ndx file - static int read_index_file(char const *filename); + int read_index_file(char const *filename); /// \brief Create atoms from a file \param filename name of the file @@ -515,13 +564,13 @@ protected: /// Output restart file colvarmodule::ofstream restart_out_os; -protected: +private: /// Counter for the current depth in the object hierarchy (useg e.g. in output) - static size_t depth_s; + size_t depth_s; /// Thread-specific depth - static std::vector depth_v; + std::vector depth_v; public: @@ -552,6 +601,10 @@ public: /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; + + /// \brief Accessor for the above + static colvarmodule *main(); + }; diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 4d2e17f098..8055d925db 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -503,7 +503,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key) std::string line; std::istringstream is(conf); - while (std::getline(is, line)) { + while (cvm::getline(is, line)) { if (line.size() == 0) continue; if (line.find_first_not_of(white_space) == @@ -539,10 +539,9 @@ int colvarparse::check_keywords(std::string &conf, char const *key) std::istream & colvarparse::getline_nocomments(std::istream &is, - std::string &line, - char const delim) + std::string &line) { - std::getline(is, line, delim); + cvm::getline(is, line); size_t const comment = line.find('#'); if (comment != std::string::npos) { line.erase(comment); diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index acdc46ddc9..9f116caafd 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -304,8 +304,7 @@ public: /// \brief Works as std::getline() but also removes everything /// between a comment character and the following newline static std::istream & getline_nocomments(std::istream &is, - std::string &s, - char const delim = '\n'); + std::string &s); /// Check if the content of the file has matching braces bool brace_check(std::string const &conf, diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index bdf65c0b1f..5b216c9d41 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -336,8 +336,6 @@ public: /// Are total forces being used? virtual bool total_forces_enabled() const { - cvm::error("Error: total forces are currently not implemented.\n", - COLVARS_NOT_IMPLEMENTED); return false; } diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index ffb83e9780..f192dcb7c0 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -12,6 +12,7 @@ #include #include "colvarscript.h" +#include "colvardeps.h" colvarscript::colvarscript(colvarproxy *p) @@ -221,6 +222,16 @@ int colvarscript::run(int argc, char const *argv[]) { } } + if (cmd == "addenergy") { + if (argc == 3) { + colvars->total_bias_energy += strtod(argv[2], NULL); + return COLVARS_OK; + } else { + result = "Wrong arguments to command \"addenergy\"\n" + help_string(); + return COLVARSCRIPT_ERROR; + } + } + result = "Syntax error\n" + help_string(); return COLVARSCRIPT_ERROR; } @@ -263,10 +274,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { } if (subcmd == "delete") { - if (cv->biases.size() > 0) { - result = "Cannot delete a colvar currently used by biases, delete those biases first"; - return COLVARSCRIPT_ERROR; + size_t i; + for (i = 0; i < cv->biases.size(); i++) { + delete cv->biases[i]; } + cv->biases.resize(0); // colvar destructor is tasked with the cleanup delete cv; // TODO this could be done by the destructors @@ -338,9 +350,8 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { return COLVARS_OK; } - if (subcmd == "state") { - cv->print_state(); - return COLVARS_OK; + if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) { + return proc_features(cv, argc, argv); } result = "Syntax error\n" + help_string(); @@ -379,11 +390,6 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { return COLVARS_OK; } - if (subcmd == "state") { - b->print_state(); - return COLVARS_OK; - } - // Subcommands for MW ABF if (subcmd == "bin") { int r = b->current_bin(); @@ -420,6 +426,10 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { return COLVARS_OK; } + if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) { + return proc_features(b, argc, argv); + } + if (argc >= 4) { std::string param = argv[3]; if (subcmd == "count") { @@ -441,6 +451,83 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { } +int colvarscript::proc_features(colvardeps *obj, + int argc, char const *argv[]) { + // size was already checked before calling + std::string subcmd = argv[2]; + + if (argc == 3) { + if (subcmd == "state") { + // TODO make this returned as result? + obj->print_state(); + return COLVARS_OK; + } + + // get and set commands require more arguments + result = "Syntax error\n" + help_string(); + return COLVARSCRIPT_ERROR; + } + + if ((subcmd == "get") || (subcmd == "set")) { + std::vector &features = obj->features(); + std::string const req_feature(argv[3]); + colvardeps::feature *f = NULL; + int fid = 0; + for (fid = 0; fid < int(features.size()); fid++) { + if (features[fid]->description == + colvarparse::to_lower_cppstr(req_feature)) { + f = features[fid]; + break; + } + } + + if (f == NULL) { + + result = "Error: feature \""+req_feature+"\" does not exist.\n"; + return COLVARSCRIPT_ERROR; + + } else { + + if (! obj->is_available(fid)) { + result = "Error: feature \""+req_feature+"\" is unavailable.\n"; + return COLVARSCRIPT_ERROR; + } + + if (subcmd == "get") { + result = cvm::to_str(obj->is_enabled(fid) ? 1 : 0); + return COLVARS_OK; + } + + if (subcmd == "set") { + if (argc == 5) { + std::string const yesno = + colvarparse::to_lower_cppstr(std::string(argv[4])); + if ((yesno == std::string("yes")) || + (yesno == std::string("on")) || + (yesno == std::string("1"))) { + obj->enable(fid); + return COLVARS_OK; + } else if ((yesno == std::string("no")) || + (yesno == std::string("off")) || + (yesno == std::string("0"))) { + // TODO disable() function does not exist yet, + // dependencies will not be resolved + // obj->disable(fid); + obj->set_enabled(fid, false); + return COLVARS_OK; + } + } + result = "Syntax error\n" + help_string(); + return COLVARSCRIPT_ERROR; + } + } + } + + result = "Syntax error\n" + help_string(); + return COLVARSCRIPT_ERROR; +} + + std::string colvarscript::help_string() { std::string buf; @@ -459,6 +546,7 @@ Input and output:\n\ load -- load a state file (requires configuration)\n\ save -- save a state file (requires configuration)\n\ update -- recalculate colvars and biases\n\ + addenergy -- add to the total bias energy\n\ printframe -- return a summary of the current frame\n\ printframelabels -- return labels to annotate printframe's output\n"; @@ -478,12 +566,17 @@ Accessing collective variables:\n\ colvar addforce -- apply given force on colvar \n\ colvar getconfig -- return config string of colvar \n\ colvar cvcflags -- enable or disable cvcs according to 0/1 flags\n\ + colvar get -- get the value of the colvar feature \n\ + colvar set -- set the value of the colvar feature \n\ \n\ Accessing biases:\n\ bias energy -- return the current energy of bias \n\ bias update -- recalculate bias \n\ bias delete -- delete bias \n\ - bias getconfig -- return config string of bias \n"; + bias getconfig -- return config string of bias \n\ + bias get -- get the value of the bias feature \n\ + bias set -- set the value of the bias feature \n\ +"; return buf; } diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index bfb8381e3b..46b1ddd203 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -51,6 +51,10 @@ private: /// Run subcommands on bias int proc_bias(int argc, char const *argv[]); + /// Run subcommands on base colvardeps object (colvar, bias, ...) + int proc_features(colvardeps *obj, + int argc, char const *argv[]); + /// Builds and return a short help std::string help_string(void); }; diff --git a/src/KOKKOS/comm_kokkos.cpp b/src/KOKKOS/comm_kokkos.cpp index 5dc1e5fa4a..f38a6c21c9 100644 --- a/src/KOKKOS/comm_kokkos.cpp +++ b/src/KOKKOS/comm_kokkos.cpp @@ -93,6 +93,9 @@ CommKokkos::~CommKokkos() void CommKokkos::init() { + maxsend = BUFMIN; + maxrecv = BUFMIN; + grow_send_kokkos(maxsend+bufextra,0,Host); grow_recv_kokkos(maxrecv,Host); diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index 7320263ba7..89453ded9b 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -352,12 +352,17 @@ int colvarproxy_lammps::smp_enabled() int colvarproxy_lammps::smp_colvars_loop() { colvarmodule *cv = this->colvars; + colvarproxy_lammps *proxy = (colvarproxy_lammps *) cv->proxy; #pragma omp parallel for - for (size_t i = 0; i < cv->colvars_smp.size(); i++) { + for (size_t i = 0; i < cv->variables_active_smp()->size(); i++) { + colvar *x = (*(cv->variables_active_smp()))[i]; + int x_item = (*(cv->variables_active_smp_items()))[i]; if (cvm::debug()) { - cvm::log("Calculating colvar \""+cv->colvars_smp[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n"); + cvm::log("["+cvm::to_str(proxy->smp_thread_id())+"/"+cvm::to_str(proxy->smp_num_threads())+ + "]: calc_colvars_items_smp(), i = "+cvm::to_str(i)+", cv = "+ + x->name+", cvc = "+cvm::to_str(x_item)+"\n"); } - cv->colvars_smp[i]->calc_cvcs(cv->colvars_smp_items[i], 1); + x->calc_cvcs(x_item, 1); } return cvm::get_error(); } @@ -367,11 +372,13 @@ int colvarproxy_lammps::smp_biases_loop() { colvarmodule *cv = this->colvars; #pragma omp parallel for - for (size_t i = 0; i < cv->biases.size(); i++) { + for (size_t i = 0; i < cv->biases_active()->size(); i++) { + colvarbias *b = (*(cv->biases_active()))[i]; if (cvm::debug()) { - cvm::log("Calculating bias \""+cv->biases[i]->name+"\" on thread "+cvm::to_str(smp_thread_id())+"\n"); + cvm::log("Calculating bias \""+b->name+"\" on thread "+ + cvm::to_str(smp_thread_id())+"\n"); } - cv->biases[i]->update(); + b->update(); } return cvm::get_error(); } diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index c0b4ea20fe..ad63eb2f9e 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -29,7 +29,7 @@ #endif #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2016-12-27" +#define COLVARPROXY_VERSION "2017-01-09" #endif /* struct for packed data communication of coordinates and forces. */ diff --git a/src/balance.cpp b/src/balance.cpp index 050f282dfe..52f6072a6c 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -456,6 +456,7 @@ void Balance::options(int iarg, int narg, char **arg) wtflag = 0; varflag = 0; + oldrcb = 0; outflag = 0; int outarg = 0; fp = NULL; @@ -491,6 +492,9 @@ void Balance::options(int iarg, int narg, char **arg) } iarg += 2+nopt; + } else if (strcmp(arg[iarg],"old") == 0) { + oldrcb = 1; + iarg++; } else if (strcmp(arg[iarg],"out") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal (fix) balance command"); outflag = 1; @@ -641,12 +645,21 @@ int *Balance::bisection(int sortflag) // invoke RCB // then invert() to create list of proc assignments for my atoms + // NOTE: (3/2017) can remove undocumented "old" option at some point + // ditto in rcb.cpp - if (wtflag) { - weight = fixstore->vstore; - rcb->compute(dim,atom->nlocal,atom->x,weight,shrinklo,shrinkhi); - } else rcb->compute(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi); - + if (oldrcb) { + if (wtflag) { + weight = fixstore->vstore; + rcb->compute_old(dim,atom->nlocal,atom->x,weight,shrinklo,shrinkhi); + } else rcb->compute_old(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi); + } else { + if (wtflag) { + weight = fixstore->vstore; + rcb->compute(dim,atom->nlocal,atom->x,weight,shrinklo,shrinkhi); + } else rcb->compute(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi); + } + rcb->invert(sortflag); // reset RCB lo/hi bounding box to full simulation box as needed diff --git a/src/balance.h b/src/balance.h index 82941b34c0..0f2f79bb15 100644 --- a/src/balance.h +++ b/src/balance.h @@ -53,6 +53,7 @@ class Balance : protected Pointers { int style; // style of LB int xflag,yflag,zflag; // xyz LB flags double *user_xsplit,*user_ysplit,*user_zsplit; // params for xyz LB + int oldrcb; // use old-style RCB compute int nitermax; // params for shift LB double stopthresh; diff --git a/src/fix_halt.cpp b/src/fix_halt.cpp index ca74efa454..2a4af4dc66 100644 --- a/src/fix_halt.cpp +++ b/src/fix_halt.cpp @@ -30,9 +30,10 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{BONDMAX,VARIABLE}; +enum{BONDMAX,TLIMIT,VARIABLE}; enum{LT,LE,GT,GE,EQ,NEQ,XOR}; enum{HARD,SOFT,CONTINUE}; +enum{NOMSG,YESMSG}; /* ---------------------------------------------------------------------- */ @@ -47,7 +48,8 @@ FixHalt::FixHalt(LAMMPS *lmp, int narg, char **arg) : idvar = NULL; - if (strcmp(arg[4],"bondmax") == 0) attribute = BONDMAX; + if (strcmp(arg[4],"tlimit") == 0) attribute = TLIMIT; + else if (strcmp(arg[4],"bondmax") == 0) attribute = BONDMAX; else if (strncmp(arg[4],"v_",2) == 0) { attribute = VARIABLE; int n = strlen(arg[4]); @@ -73,6 +75,7 @@ FixHalt::FixHalt(LAMMPS *lmp, int narg, char **arg) : // parse optional args eflag = SOFT; + msgflag = YESMSG; int iarg = 7; while (iarg < narg) { @@ -83,6 +86,12 @@ FixHalt::FixHalt(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg+1],"continue") == 0) eflag = CONTINUE; else error->all(FLERR,"Illegal fix halt command"); iarg += 2; + } else if (strcmp(arg[iarg],"message") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix halt command"); + if (strcmp(arg[iarg+1],"no") == 0) msgflag = NOMSG; + else if (strcmp(arg[iarg+1],"yes") == 0) msgflag = YESMSG; + else error->all(FLERR,"Illegal fix halt command"); + iarg += 2; } else error->all(FLERR,"Illegal fix halt command"); } @@ -125,6 +134,11 @@ void FixHalt::init() if (input->variable->equalstyle(ivar) == 0) error->all(FLERR,"Fix halt variable is not equal-style variable"); } + + // settings used by TLIMIT + + nextstep = (update->ntimestep/nevery)*nevery + nevery; + tratio = 0.5; } /* ---------------------------------------------------------------------- */ @@ -135,8 +149,12 @@ void FixHalt::end_of_step() double attvalue; - if (attribute == BONDMAX) attvalue = bondmax(); - else { + if (attribute == TLIMIT) { + if (update->ntimestep != nextstep) return; + attvalue = tlimit(); + } else if (attribute == BONDMAX) { + attvalue = bondmax(); + } else { modify->clearstep_compute(); attvalue = input->variable->compute_equal(ivar); modify->addstep_compute(update->ntimestep + nevery); @@ -169,9 +187,10 @@ void FixHalt::end_of_step() sprintf(str,"Fix halt %s condition met on step %ld with value %g", id,update->ntimestep,attvalue); - if (eflag == HARD) error->all(FLERR,str); - else if (eflag == SOFT || eflag == CONTINUE) { - if (comm->me == 0) error->message(FLERR,str); + if (eflag == HARD) { + error->all(FLERR,str); + } else if (eflag == SOFT || eflag == CONTINUE) { + if (comm->me == 0 && msgflag == YESMSG) error->message(FLERR,str); timer->force_timeout(); } } @@ -218,3 +237,27 @@ double FixHalt::bondmax() return sqrt(maxall); } + +/* ---------------------------------------------------------------------- + compute synced elapsed time + reset nextstep = estimate of timestep when run will end + first project to 1/2 the run time, thereafter to end of run +------------------------------------------------------------------------- */ + +double FixHalt::tlimit() +{ + double cpu = timer->elapsed(Timer::TOTAL); + MPI_Bcast(&cpu,1,MPI_DOUBLE,0,world); + + if (cpu < value) { + bigint elapsed = update->ntimestep - update->firststep; + bigint final = update->firststep + + static_cast (tratio*value/cpu * elapsed); + nextstep = (final/nevery)*nevery + nevery; + tratio = 1.0; + } + + //printf("EVAL %ld %g %d\n",update->ntimestep,cpu,nevery); + + return cpu; +} diff --git a/src/fix_halt.h b/src/fix_halt.h index 156397a4b5..52d2f95d2f 100644 --- a/src/fix_halt.h +++ b/src/fix_halt.h @@ -35,11 +35,13 @@ class FixHalt : public Fix { void post_run(); private: - int attribute,operation,eflag,ivar; - double value; + int attribute,operation,eflag,msgflag,ivar; + bigint nextstep; + double value,tratio; char *idvar; double bondmax(); + double tlimit(); }; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 48182d9f0b..759469eae4 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1129,6 +1129,11 @@ void Neighbor::morph_halffull() if (!irq->half) continue; + // Kokkos doesn't yet support half from full + + if (irq->kokkos_host) continue; + if (irq->kokkos_device) continue; + // these lists are created other ways, no need for halffull // do want to process skip lists @@ -1154,8 +1159,6 @@ void Neighbor::morph_halffull() // this includes custom cutoff set by requestor // no need to check respaouter b/c it stores same pairs // no need to check dnum b/c only set for history - // NOTE: need check for 2 Kokkos flags? - // Kokkos doesn't yet support half from full? if (irq->ghost != jrq->ghost) continue; if (irq->size != jrq->size) continue; diff --git a/src/rcb.cpp b/src/rcb.cpp index 7223e4616a..1aca015daf 100644 --- a/src/rcb.cpp +++ b/src/rcb.cpp @@ -42,7 +42,7 @@ RCB::RCB(LAMMPS *lmp) : Pointers(lmp) dots = NULL; nlist = maxlist = 0; - dotlist = dotmark = NULL; + dotlist = dotmark = dotmark_select = NULL; maxbuf = 0; buf = NULL; @@ -73,6 +73,7 @@ RCB::~RCB() memory->sfree(dots); memory->destroy(dotlist); memory->destroy(dotmark); + memory->destroy(dotmark_select); memory->sfree(buf); memory->destroy(recvproc); @@ -91,6 +92,9 @@ RCB::~RCB() /* ---------------------------------------------------------------------- perform RCB balancing of N particles at coords X in bounding box LO/HI + NEW version: each RCB cut is tested in all dimensions + dimeension that produces 2 boxes with largest min size is selected + this is to prevent very narrow boxes from being produced if wt = NULL, ignore per-particle weights if wt defined, per-particle weights > 0.0 dimension = 2 or 3 @@ -103,6 +107,523 @@ RCB::~RCB() void RCB::compute(int dimension, int n, double **x, double *wt, double *bboxlo, double *bboxhi) +{ + int i,j,k; + int keep,outgoing,incoming,incoming2; + int dim,markactive; + int indexlo,indexhi; + int first_iteration,breakflag; + double wttot,wtlo,wthi,wtsum,wtok,wtupto,wtmax; + double targetlo,targethi; + double valuemin,valuemax,valuehalf,valuehalf_select,smaller; + double tolerance; + MPI_Comm comm,comm_half; + MPI_Request request,request2; + Median med,medme; + + // create list of my Dots + + ndot = nkeep = noriginal = n; + + if (ndot > maxdot) { + maxdot = ndot; + memory->sfree(dots); + dots = (Dot *) memory->smalloc(ndot*sizeof(Dot),"RCB:dots"); + } + + for (i = 0; i < ndot; i++) { + dots[i].x[0] = x[i][0]; + dots[i].x[1] = x[i][1]; + dots[i].x[2] = x[i][2]; + dots[i].proc = me; + dots[i].index = i; + } + + if (wt) + for (i = 0; i < ndot; i++) dots[i].wt = wt[i]; + else + for (i = 0; i < ndot; i++) dots[i].wt = 1.0; + + // initial bounding box = simulation box + // includes periodic or shrink-wrapped boundaries + + lo = bbox.lo; + hi = bbox.hi; + + lo[0] = bboxlo[0]; + lo[1] = bboxlo[1]; + lo[2] = bboxlo[2]; + hi[0] = bboxhi[0]; + hi[1] = bboxhi[1]; + hi[2] = bboxhi[2]; + + cut = 0.0; + cutdim = -1; + + // initialize counters + + counters[0] = 0; + counters[1] = 0; + counters[2] = 0; + counters[3] = ndot; + counters[4] = maxdot; + counters[5] = 0; + counters[6] = 0; + + // create communicator for use in recursion + + MPI_Comm_dup(world,&comm); + + // recurse until partition is a single proc = me + // proclower,procupper = lower,upper procs in partition + // procmid = 1st proc in upper half of partition + + int procpartner,procpartner2; + + int procmid; + int proclower = 0; + int procupper = nprocs - 1; + + while (proclower != procupper) { + + // if odd # of procs, lower partition gets extra one + + procmid = proclower + (procupper - proclower) / 2 + 1; + + // determine communication partner(s) + // readnumber = # of proc partners to read from + + if (me < procmid) + procpartner = me + (procmid - proclower); + else + procpartner = me - (procmid - proclower); + + int readnumber = 1; + if (procpartner > procupper) { + readnumber = 0; + procpartner--; + } + if (me == procupper && procpartner != procmid - 1) { + readnumber = 2; + procpartner2 = procpartner + 1; + } + + // wttot = summed weight of entire partition + // search tolerance = largest single weight (plus epsilon) + // targetlo = desired weight in lower half of partition + // targethi = desired weight in upper half of partition + + wtmax = wtsum = 0.0; + + if (wt) { + for (i = 0; i < ndot; i++) { + wtsum += dots[i].wt; + if (dots[i].wt > wtmax) wtmax = dots[i].wt; + } + } else { + for (i = 0; i < ndot; i++) wtsum += dots[i].wt; + wtmax = 1.0; + } + + MPI_Allreduce(&wtsum,&wttot,1,MPI_DOUBLE,MPI_SUM,comm); + if (wt) MPI_Allreduce(&wtmax,&tolerance,1,MPI_DOUBLE,MPI_MAX,comm); + else tolerance = 1.0; + + tolerance *= 1.0 + TINY; + targetlo = wttot * (procmid - proclower) / (procupper + 1 - proclower); + targethi = wttot - targetlo; + + // attempt a cut in each dimension + // each cut produces 2 boxes, each with a reduced box length in that dim + // smaller = the smaller of the 2 reduced box lengths in that dimension + // choose to cut in dimension which produces largest smaller value + // this should induce final proc sub-boxes to be as cube-ish as possible + // dim_select = selected cut dimension + // valuehalf_select = valuehalf in that dimension + // dotmark_select = dot markings in that dimension + + int dim_select = -1; + double largest = 0.0; + + for (dim = 0; dim < dimension; dim++) { + + // create active list and mark array for dots + // initialize active list to all dots + + if (ndot > maxlist) { + memory->destroy(dotlist); + memory->destroy(dotmark); + memory->destroy(dotmark_select); + maxlist = maxdot; + memory->create(dotlist,maxlist,"RCB:dotlist"); + memory->create(dotmark,maxlist,"RCB:dotmark"); + memory->create(dotmark_select,maxlist,"RCB:dotmark_select"); + } + + nlist = ndot; + for (i = 0; i < nlist; i++) dotlist[i] = i; + + // median iteration + // zoom in on bisector until correct # of dots in each half of partition + // as each iteration of median-loop begins, require: + // all non-active dots are marked with 0/1 in dotmark + // valuemin <= every active dot <= valuemax + // wtlo, wthi = total wt of non-active dots + // when leave median-loop, require only: + // valuehalf = correct cut position + // all dots <= valuehalf are marked with 0 in dotmark + // all dots >= valuehalf are marked with 1 in dotmark + // markactive = which side of cut is active = 0/1 + // indexlo,indexhi = indices of dot closest to median + + wtlo = wthi = 0.0; + valuemin = lo[dim]; + valuemax = hi[dim]; + first_iteration = 1; + indexlo = indexhi = 0; + + while (1) { + + // choose bisector value + // use old value on 1st iteration if old cut dimension is the same + // on 2nd option: could push valuehalf towards geometric center + // with "1.0-factor" to force overshoot + + if (first_iteration && reuse && dim == tree[procmid].dim) { + counters[5]++; + valuehalf = tree[procmid].cut; + if (valuehalf < valuemin || valuehalf > valuemax) + valuehalf = 0.5 * (valuemin + valuemax); + } else if (wt) + valuehalf = valuemin + (targetlo - wtlo) / + (wttot - wtlo - wthi) * (valuemax - valuemin); + else + valuehalf = 0.5 * (valuemin + valuemax); + + first_iteration = 0; + + // initialize local median data structure + + medme.totallo = medme.totalhi = 0.0; + medme.valuelo = -MYHUGE; + medme.valuehi = MYHUGE; + medme.wtlo = medme.wthi = 0.0; + medme.countlo = medme.counthi = 0; + medme.proclo = medme.prochi = me; + + // mark all active dots on one side or other of bisector + // also set all fields in median data struct + // save indices of closest dots on either side + + for (j = 0; j < nlist; j++) { + i = dotlist[j]; + if (dots[i].x[dim] <= valuehalf) { // in lower part + medme.totallo += dots[i].wt; + dotmark[i] = 0; + if (dots[i].x[dim] > medme.valuelo) { // my closest dot + medme.valuelo = dots[i].x[dim]; + medme.wtlo = dots[i].wt; + medme.countlo = 1; + indexlo = i; + } else if (dots[i].x[dim] == medme.valuelo) { // tied for closest + medme.wtlo += dots[i].wt; + medme.countlo++; + } + } + else { // in upper part + medme.totalhi += dots[i].wt; + dotmark[i] = 1; + if (dots[i].x[dim] < medme.valuehi) { // my closest dot + medme.valuehi = dots[i].x[dim]; + medme.wthi = dots[i].wt; + medme.counthi = 1; + indexhi = i; + } else if (dots[i].x[dim] == medme.valuehi) { // tied for closest + medme.wthi += dots[i].wt; + medme.counthi++; + } + } + } + + // combine median data struct across current subset of procs + + counters[0]++; + MPI_Allreduce(&medme,&med,1,med_type,med_op,comm); + + // test median guess for convergence + // move additional dots that are next to cut across it + + if (wtlo + med.totallo < targetlo) { // lower half TOO SMALL + + wtlo += med.totallo; + valuehalf = med.valuehi; + + if (med.counthi == 1) { // only one dot to move + if (wtlo + med.wthi < targetlo) { // move it, keep iterating + if (me == med.prochi) dotmark[indexhi] = 0; + } + else { // only move if beneficial + if (wtlo + med.wthi - targetlo < targetlo - wtlo) + if (me == med.prochi) dotmark[indexhi] = 0; + break; // all done + } + } + else { // multiple dots to move + breakflag = 0; + wtok = 0.0; + if (medme.valuehi == med.valuehi) wtok = medme.wthi; + if (wtlo + med.wthi >= targetlo) { // all done + MPI_Scan(&wtok,&wtupto,1,MPI_DOUBLE,MPI_SUM,comm); + wtmax = targetlo - wtlo; + if (wtupto > wtmax) wtok = wtok - (wtupto - wtmax); + breakflag = 1; + } // wtok = most I can move + for (j = 0, wtsum = 0.0; j < nlist && wtsum < wtok; j++) { + i = dotlist[j]; + if (dots[i].x[dim] == med.valuehi) { // only move if better + if (wtsum + dots[i].wt - wtok < wtok - wtsum) + dotmark[i] = 0; + wtsum += dots[i].wt; + } + } + if (breakflag) break; // done if moved enough + } + + wtlo += med.wthi; + if (targetlo-wtlo <= tolerance) break; // close enough + + valuemin = med.valuehi; // iterate again + markactive = 1; + } + + else if (wthi + med.totalhi < targethi) { // upper half TOO SMALL + + wthi += med.totalhi; + valuehalf = med.valuelo; + + if (med.countlo == 1) { // only one dot to move + if (wthi + med.wtlo < targethi) { // move it, keep iterating + if (me == med.proclo) dotmark[indexlo] = 1; + } + else { // only move if beneficial + if (wthi + med.wtlo - targethi < targethi - wthi) + if (me == med.proclo) dotmark[indexlo] = 1; + break; // all done + } + } + else { // multiple dots to move + breakflag = 0; + wtok = 0.0; + if (medme.valuelo == med.valuelo) wtok = medme.wtlo; + if (wthi + med.wtlo >= targethi) { // all done + MPI_Scan(&wtok,&wtupto,1,MPI_DOUBLE,MPI_SUM,comm); + wtmax = targethi - wthi; + if (wtupto > wtmax) wtok = wtok - (wtupto - wtmax); + breakflag = 1; + } // wtok = most I can move + for (j = 0, wtsum = 0.0; j < nlist && wtsum < wtok; j++) { + i = dotlist[j]; + if (dots[i].x[dim] == med.valuelo) { // only move if better + if (wtsum + dots[i].wt - wtok < wtok - wtsum) + dotmark[i] = 1; + wtsum += dots[i].wt; + } + } + if (breakflag) break; // done if moved enough + } + + wthi += med.wtlo; + if (targethi-wthi <= tolerance) break; // close enough + + valuemax = med.valuelo; // iterate again + markactive = 0; + } + + else // Goldilocks result: both partitions just right + break; + + // shrink the active list + + k = 0; + for (j = 0; j < nlist; j++) { + i = dotlist[j]; + if (dotmark[i] == markactive) dotlist[k++] = i; + } + nlist = k; + } + + // cut produces 2 sub-boxes with reduced size in dim + // compare smaller of the 2 sizes to previous dims + // keep dim that has the largest smaller + + smaller = MIN(valuehalf-lo[dim],hi[dim]-valuehalf); + if (smaller > largest) { + largest = smaller; + dim_select = dim; + valuehalf_select = valuehalf; + memcpy(dotmark_select,dotmark,ndot*sizeof(int)); + } + } + + // copy results for best dim cut into dim,valuehalf,dotmark + + dim = dim_select; + valuehalf = valuehalf_select; + memcpy(dotmark,dotmark_select,ndot*sizeof(int)); + + // found median + // store cut info only if I am procmid + + if (me == procmid) { + cut = valuehalf; + cutdim = dim; + } + + // use cut to shrink my RCB bounding box + + if (me < procmid) hi[dim] = valuehalf; + else lo[dim] = valuehalf; + + // outgoing = number of dots to ship to partner + // nkeep = number of dots that have never migrated + + markactive = (me < procpartner); + for (i = 0, keep = 0, outgoing = 0; i < ndot; i++) + if (dotmark[i] == markactive) outgoing++; + else if (i < nkeep) keep++; + nkeep = keep; + + // alert partner how many dots I'll send, read how many I'll recv + + MPI_Send(&outgoing,1,MPI_INT,procpartner,0,world); + incoming = 0; + if (readnumber) { + MPI_Recv(&incoming,1,MPI_INT,procpartner,0,world,MPI_STATUS_IGNORE); + if (readnumber == 2) { + MPI_Recv(&incoming2,1,MPI_INT,procpartner2,0,world,MPI_STATUS_IGNORE); + incoming += incoming2; + } + } + + // check if need to alloc more space + + int ndotnew = ndot - outgoing + incoming; + if (ndotnew > maxdot) { + while (maxdot < ndotnew) maxdot += DELTA; + dots = (Dot *) memory->srealloc(dots,maxdot*sizeof(Dot),"RCB::dots"); + counters[6]++; + } + + counters[1] += outgoing; + counters[2] += incoming; + if (ndotnew > counters[3]) counters[3] = ndotnew; + if (maxdot > counters[4]) counters[4] = maxdot; + + // malloc comm send buffer + + if (outgoing > maxbuf) { + memory->sfree(buf); + maxbuf = outgoing; + buf = (Dot *) memory->smalloc(maxbuf*sizeof(Dot),"RCB:buf"); + } + + // fill buffer with dots that are marked for sending + // pack down the unmarked ones + + keep = outgoing = 0; + for (i = 0; i < ndot; i++) { + if (dotmark[i] == markactive) + memcpy(&buf[outgoing++],&dots[i],sizeof(Dot)); + else + memcpy(&dots[keep++],&dots[i],sizeof(Dot)); + } + + // post receives for dots + + if (readnumber > 0) { + MPI_Irecv(&dots[keep],incoming*sizeof(Dot),MPI_CHAR, + procpartner,1,world,&request); + if (readnumber == 2) { + keep += incoming - incoming2; + MPI_Irecv(&dots[keep],incoming2*sizeof(Dot),MPI_CHAR, + procpartner2,1,world,&request2); + } + } + + // handshake before sending dots to insure recvs have been posted + + if (readnumber > 0) { + MPI_Send(NULL,0,MPI_INT,procpartner,0,world); + if (readnumber == 2) MPI_Send(NULL,0,MPI_INT,procpartner2,0,world); + } + MPI_Recv(NULL,0,MPI_INT,procpartner,0,world,MPI_STATUS_IGNORE); + + // send dots to partner + + MPI_Rsend(buf,outgoing*sizeof(Dot),MPI_CHAR,procpartner,1,world); + + // wait until all dots are received + + if (readnumber > 0) { + MPI_Wait(&request,MPI_STATUS_IGNORE); + if (readnumber == 2) MPI_Wait(&request2,MPI_STATUS_IGNORE); + } + + ndot = ndotnew; + + // cut partition in half, create new communicators of 1/2 size + + int split; + if (me < procmid) { + procupper = procmid - 1; + split = 0; + } else { + proclower = procmid; + split = 1; + } + + MPI_Comm_split(comm,split,me,&comm_half); + MPI_Comm_free(&comm); + comm = comm_half; + } + + // clean up + + MPI_Comm_free(&comm); + + // set public variables with results of rebalance + + nfinal = ndot; + + if (nfinal > maxrecv) { + memory->destroy(recvproc); + memory->destroy(recvindex); + maxrecv = nfinal; + memory->create(recvproc,maxrecv,"RCB:recvproc"); + memory->create(recvindex,maxrecv,"RCB:recvindex"); + } + + for (i = 0; i < nfinal; i++) { + recvproc[i] = dots[i].proc; + recvindex[i] = dots[i].index; + } +} + +/* ---------------------------------------------------------------------- + perform RCB balancing of N particles at coords X in bounding box LO/HI + OLD version: each RCB cut is made in longest dimension of sub-box + if wt = NULL, ignore per-particle weights + if wt defined, per-particle weights > 0.0 + dimension = 2 or 3 + as documented in rcb.h: + sets noriginal,nfinal,nkeep,recvproc,recvindex,lo,hi + all proc particles will be inside or on surface of 3-d box + defined by final lo/hi + // NOTE: worry about re-use of data structs for fix balance? +------------------------------------------------------------------------- */ + +void RCB::compute_old(int dimension, int n, double **x, double *wt, + double *bboxlo, double *bboxhi) { int i,j,k; int keep,outgoing,incoming,incoming2; diff --git a/src/rcb.h b/src/rcb.h index 7a82d42f93..90b82f5952 100644 --- a/src/rcb.h +++ b/src/rcb.h @@ -42,6 +42,7 @@ class RCB : protected Pointers { RCB(class LAMMPS *); ~RCB(); void compute(int, int, double **, double *, double *, double *); + void compute_old(int, int, double **, double *, double *, double *); void invert(int sortflag = 0); bigint memory_usage(); @@ -99,6 +100,7 @@ class RCB : protected Pointers { int maxlist; int *dotlist; int *dotmark; + int *dotmark_select; int maxbuf; Dot *buf; diff --git a/src/thermo.cpp b/src/thermo.cpp index 75e72ada64..dbbeff4998 100644 --- a/src/thermo.cpp +++ b/src/thermo.cpp @@ -55,7 +55,8 @@ using namespace MathConst; // customize a new keyword by adding to this list: -// step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, part, timeremain +// step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, +// part, timeremain // atoms, temp, press, pe, ke, etotal, enthalpy // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail // vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, @@ -394,6 +395,11 @@ void Thermo::compute(int flag) if (flushflag) fflush(logfile); } } + + // set to 1, so that subsequent invocations of CPU time will be non-zero + // e.g. via variables in print command + + firststep = 1; } /* ---------------------------------------------------------------------- diff --git a/src/version.h b/src/version.h index 703db28c26..e0aa18371c 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "7 Mar 2017" +#define LAMMPS_VERSION "10 Mar 2017"