ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into USER-DPD_kokkos as of patch 10 Mar 2017

Browse Source
This commit is contained in:
Tim Mattox
2017-03-10 21:37:52 -05:00
parent c7eadfdefd f871ecdc67
commit 17c8e695a4
45 changed files with 2063 additions and 919 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/Manual.txt
View File

@ -1,7 +1,7 @@
<!-- HTML_ONLY --> <!-- HTML_ONLY -->
<HEAD> <HEAD>
<TITLE>LAMMPS Users Manual</TITLE> <TITLE>LAMMPS Users Manual</TITLE>
<META NAME="docnumber" CONTENT="7 Mar 2017 version"> <META NAME="docnumber" CONTENT="10 Mar 2017 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories"> <META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License."> <META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD> </HEAD>
@ -21,7 +21,7 @@
<H1></H1> <H1></H1>
LAMMPS Documentation :c,h3 LAMMPS Documentation :c,h3
7 Mar 2017 version :c,h4 10 Mar 2017 version :c,h4
Version info: :h4 Version info: :h4

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

38
doc/src/balance.txt
View File

@ -286,24 +286,32 @@ above. It performs a recursive coordinate bisectioning (RCB) of the
simulation domain. The basic idea is as follows. simulation domain. The basic idea is as follows.
The simulation domain is cut into 2 boxes by an axis-aligned cut in 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. one of the dimensions, leaving one new sub-box on either side of the
All the processors are also partitioned into 2 groups, half assigned cut. Which dimension is chosen for the cut depends on the particle
to the box on the lower side of the cut, and half to the box on the (weight) distribution within the parent box. Normally the longest
upper side. (If the processor count is odd, one side gets an extra dimension of the box is cut, but if all (or most) of the particles are
processor.) The cut is positioned so that the number of particles in at one end of the box, a cut may be performed in another dimension to
the lower box is exactly the number that the processors assigned to induce sub-boxes that are more cube-ish (3d) or square-ish (2d) in
that box should own for load balance to be perfect. This also makes shape.
load balance for the upper box perfect. The positioning is done
iteratively, by a bisectioning method. Note that counting particles After the cut is made, all the processors are also partitioned into 2
on either side of the cut requires communication between all groups, half assigned to the box on the lower side of the cut, and
processors at each iteration. 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 That is the procedure for the first cut. Subsequent cuts are made
recursively, in exactly the same manner. The subset of processors recursively, in exactly the same manner. The subset of processors
assigned to each box make a new cut in the longest dimension of that assigned to each box make a new cut in one dimension of that box,
box, splitting the box, the subset of processors, and the particles splitting the box, the subset of processors, and the particles in the
in the box in two. The recursion continues until every processor is box in two. The recursion continues until every processor is assigned
assigned a sub-box of the entire simulation domain, and owns the a sub-box of the entire simulation domain, and owns the (weighted)
particles in that sub-box. particles in that sub-box.
:line :line

4
doc/src/change_box.txt
View File

@ -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: is non-periodic and atoms extend from 0 to 20 in all dimensions:
change_box all x final -10 20 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 change_box all x final -10 20 boundary f s s
create_atoms 1 single -5 5 5 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 NOTE: Unlike the earlier "displace_box" version of this command, atom
remapping is NOT performed by default. This command allows remapping remapping is NOT performed by default. This command allows remapping

11
doc/src/create_atoms.txt
View File

@ -134,6 +134,17 @@ not overlap existing atoms inappropriately, especially if molecules
are being added. The "delete_atoms"_delete_atoms.html command can be are being added. The "delete_atoms"_delete_atoms.html command can be
used to remove overlapping atoms or molecules. 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 :line
Individual atoms are inserted by this command, unless the {mol} Individual atoms are inserted by this command, unless the {mol}

54
doc/src/fix_halt.txt
View File

@ -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 ID, group-ID are documented in "fix"_fix.html command :ulb,l
halt = style name of this fix command :l halt = style name of this fix command :l
N = check halt condition every N steps :l N = check halt condition every N steps :l
attribute = hstyle or v_name :l attribute = {bondmax} or {tlimit} or v_name :l
hstyle = {bondmax} bondmax = length of longest bond in the system
tlimit = elapsed CPU time
v_name = name of "equal-style variable"_variable.html :pre v_name = name of "equal-style variable"_variable.html :pre
operator = "<" or "<=" or ">" or ">=" or "==" or "!=" or "|^" :l operator = "<" or "<=" or ">" or ">=" or "==" or "!=" or "|^" :l
avalue = numeric value to compare attribute to :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 zero or more keyword/value pairs may be appended :l
keyword = {error} :l keyword = {error} or {message} :l
{error} value = {hard} or {soft} or {continue} :pre {error} value = {hard} or {soft} or {continue}
{message} value = {yes} or {no} :pre
:ule :ule
[Examples:] [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 group-ID is ignored by this fix.
The specified {attribute} can be one of the {hstyle} options listed The specified {attribute} can be one of the options listed above,
above, or an "equal-style variable"_variable.html referenced as namely {bondmax} or {tlimit}, or an "equal-style
{v_name}, where "name" is the name of a variable that has been defined variable"_variable.html referenced as {v_name}, where "name" is the
previously in the input script. name of a variable that has been defined previously in the input
script.
The only {hstyle} option currently implemented is {bondmax}. This The {bondmax} attribute will loop over all bonds in the system,
will loop over all bonds in the system, compute their current compute their current lengths, and set {attribute} to the longest bond
lengths, and set {attribute} to the longest bond distance. 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 Equal-style variables evaluate to a numeric value. See the
"variable"_variable.html command for a description. They calculate "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 halt ID, so that the same condition is not immediately triggered in a
subsequent run. 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:] [Restart, fix_modify, output, run start/stop, minimize info:]
No information about this fix is written to "binary restart 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:] [Default:]
The option defaults are error = hard. The option defaults are error = hard and message = yes.

2
doc/src/quit.txt
View File

@ -17,7 +17,7 @@ status = numerical exit status (optional)
[Examples:] [Examples:]
quit quit
if "$n > 10000" then "quit 1":pre if "$n > 10000" then "quit 1" :pre
[Description:] [Description:]

2
doc/src/thermo_style.txt
View File

@ -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 elaplong = timesteps since start of initial run in a series of runs
dt = timestep size dt = timestep size
time = simulation time 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 tpcpu = time per CPU second
spcpu = timesteps per CPU second spcpu = timesteps per CPU second
cpuremain = estimated CPU time remaining in run cpuremain = estimated CPU time remaining in run

362
lib/colvars/colvar.cpp
View File

@ -1,5 +1,5 @@
// -*- c++ -*-
// -*- c++ -*-
// This file is part of the Collective Variables module (Colvars). // This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at: // The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars // https://github.com/colvars/colvars
@ -7,6 +7,7 @@
// If you wish to distribute your changes, please submit them to the // If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub. // Colvars repository at GitHub.
#include "colvarmodule.h" #include "colvarmodule.h"
#include "colvarvalue.h" #include "colvarvalue.h"
#include "colvarparse.h" #include "colvarparse.h"
@ -23,20 +24,31 @@ bool compare(colvar::cvc *i, colvar::cvc *j) {
} }
colvar::colvar(std::string const &conf) colvar::colvar()
: colvarparse(conf) {
// 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"); cvm::log("Initializing a new collective variable.\n");
colvarparse::init(conf);
int error_code = COLVARS_OK; int error_code = COLVARS_OK;
get_keyval(conf, "name", this->name, colvarmodule *cv = cvm::main();
(std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));
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+ cvm::error("Error: this colvar cannot have the same name, \""+this->name+
"\", as another colvar.\n", "\", as another colvar.\n",
INPUT_ERROR); INPUT_ERROR);
return; return INPUT_ERROR;
} }
// Initialize dependency members // Initialize dependency members
@ -44,14 +56,13 @@ colvar::colvar(std::string const &conf)
this->description = "colvar " + this->name; this->description = "colvar " + this->name;
// Initialize static array once and for all
init_cv_requires();
kinetic_energy = 0.0; kinetic_energy = 0.0;
potential_energy = 0.0; potential_energy = 0.0;
error_code |= init_components(conf); error_code |= init_components(conf);
if (error_code != COLVARS_OK) return; if (error_code != COLVARS_OK) {
return cvm::get_error();
}
size_t i; size_t i;
@ -59,8 +70,6 @@ colvar::colvar(std::string const &conf)
if (get_keyval(conf, "scriptedFunction", scripted_function, if (get_keyval(conf, "scriptedFunction", scripted_function,
"", colvarparse::parse_silent)) { "", colvarparse::parse_silent)) {
// Make feature available only on user request
provide(f_cv_scripted);
enable(f_cv_scripted); enable(f_cv_scripted);
cvm::log("This colvar uses scripted function \"" + scripted_function + "\"."); 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) { if (x.type() == colvarvalue::type_notset) {
cvm::error("Could not parse scripted colvar type."); cvm::error("Could not parse scripted colvar type.", INPUT_ERROR);
return; return INPUT_ERROR;
} }
cvm::log(std::string("Expecting colvar value of type ") 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) { if (x.type() == colvarvalue::type_vector) {
int size; int size;
if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) { if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
cvm::error("Error: no size specified for vector scripted function."); cvm::error("Error: no size specified for vector scripted function.",
return; INPUT_ERROR);
return INPUT_ERROR;
} }
x.vector1d_value.resize(size); 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: // Colvar is homogeneous if:
@ -168,7 +178,7 @@ colvar::colvar(std::string const &conf)
homogeneous = false; homogeneous = false;
} }
} }
feature_states[f_cv_homogeneous]->enabled = homogeneous; set_enabled(f_cv_homogeneous, homogeneous);
} }
// Colvar is deemed periodic if: // Colvar is deemed periodic if:
@ -188,7 +198,7 @@ colvar::colvar(std::string const &conf)
"Make sure that you know what you are doing!"); "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 // check that cvcs are compatible
@ -202,7 +212,7 @@ colvar::colvar(std::string const &conf)
"by using components of different types. " "by using components of different types. "
"You must use the same type in order to " "You must use the same type in order to "
"sum them together.\n", INPUT_ERROR); "sum them together.\n", INPUT_ERROR);
return; return INPUT_ERROR;
} }
} }
@ -215,44 +225,110 @@ colvar::colvar(std::string const &conf)
f.type(value()); f.type(value());
f_accumulated.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(); 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); get_keyval(conf, "width", width, 1.0);
if (width <= 0.0) { if (width <= 0.0) {
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
return; return INPUT_ERROR;
} }
lower_boundary.type(value()); lower_boundary.type(value());
lower_wall.type(value());
upper_boundary.type(value()); upper_boundary.type(value());
upper_wall.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 (is_enabled(f_cv_scalar)) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) { if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
provide(f_cv_lower_boundary);
enable(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 (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) {
if (lower_wall_k > 0.0) { 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); 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)) { if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
provide(f_cv_upper_boundary);
enable(f_cv_upper_boundary); enable(f_cv_upper_boundary);
} }
get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0); if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) {
if (upper_wall_k > 0.0) { 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); 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, "+ ", is not higher than the lower boundary, "+
cvm::to_str(lower_boundary)+".\n", cvm::to_str(lower_boundary)+".\n",
INPUT_ERROR); INPUT_ERROR);
} return 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);
} }
} }
@ -289,83 +356,90 @@ colvar::colvar(std::string const &conf)
cvm::error("Error: trying to expand boundaries that already " cvm::error("Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.\n", "cover a whole period of a periodic colvar.\n",
INPUT_ERROR); INPUT_ERROR);
return INPUT_ERROR;
} }
if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) { if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
cvm::error("Error: inconsistent configuration " cvm::error("Error: inconsistent configuration "
"(trying to expand boundaries with both " "(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).\n", "hardLowerBoundary and hardUpperBoundary enabled).\n",
INPUT_ERROR); 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) { int colvar::init_extended_Lagrangian(std::string const &conf)
cvm::real temp, tolerance, period; {
bool b_extended_Lagrangian;
get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
cvm::log("Enabling the extended Lagrangian term for colvar \""+ if (b_extended_Lagrangian) {
this->name+"\".\n"); cvm::real temp, tolerance, period;
// Make feature available only on user request cvm::log("Enabling the extended Lagrangian term for colvar \""+
provide(f_cv_extended_Lagrangian); this->name+"\".\n");
enable(f_cv_extended_Lagrangian);
provide(f_cv_Langevin);
// The extended mass will apply forces enable(f_cv_extended_Lagrangian);
enable(f_cv_gradient);
xr.type(value()); xr.type(value());
vr.type(value()); vr.type(value());
fr.type(value()); fr.type(value());
const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature());
if (temp <= 0.0) { if (temp <= 0.0) {
if (found) if (found)
cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR); cvm::error("Error: \"extendedTemp\" must be positive.\n", INPUT_ERROR);
else else
cvm::error("Error: a positive temperature must be provided, either " cvm::error("Error: a positive temperature must be provided, either "
"by enabling a thermostat, or through \"extendedTemp\".\n", "by enabling a thermostat, or through \"extendedTemp\".\n",
INPUT_ERROR); 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); get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
if (tolerance <= 0.0) { if (ext_gamma < 0.0) {
cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR); cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\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"); if (ext_gamma != 0.0) {
enable(f_cv_Langevin);
get_keyval(conf, "extendedTimeConstant", period, 200.0); ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1
if (period <= 0.0) { ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());
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());
}
} }
} }
return COLVARS_OK;
}
int colvar::init_output_flags(std::string const &conf)
{
{ {
bool b_output_value; bool b_output_value;
get_keyval(conf, "outputValue", b_output_value, true); 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)) { if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n" 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); "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, "outputAppliedForce", f_cv_output_applied_force, false);
get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false); get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
// Start in active state by default return COLVARS_OK;
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");
} }
@ -637,7 +690,7 @@ int colvar::parse_analysis(std::string const &conf)
std::string runave_outfile; std::string runave_outfile;
get_keyval(conf, "runAveOutputFile", runave_outfile, get_keyval(conf, "runAveOutputFile", runave_outfile,
std::string(cvm::output_prefix+"."+ std::string(cvm::output_prefix()+"."+
this->name+".runave.traj")); this->name+".runave.traj"));
size_t const this_cv_width = x.output_width(cvm::cv_width); 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, "corrFuncNormalize", acf_normalize, true);
get_keyval(conf, "corrFuncOutputFile", acf_outfile, get_keyval(conf, "corrFuncOutputFile", acf_outfile,
std::string(cvm::output_prefix+"."+this->name+ std::string(cvm::output_prefix()+"."+this->name+
".corrfunc.dat")); ".corrfunc.dat"));
} }
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -730,11 +783,12 @@ colvar::~colvar()
} }
// remove reference to this colvar from the CVM // remove reference to this colvar from the CVM
for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin(); colvarmodule *cv = cvm::main();
cvi != cvm::colvars.end(); for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
cvi != cv->variables()->end();
++cvi) { ++cvi) {
if ( *cvi == this) { if ( *cvi == this) {
cvm::colvars.erase(cvi); cv->variables()->erase(cvi);
break; break;
} }
} }
@ -892,7 +946,6 @@ int colvar::collect_cvc_values()
cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n"); cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
if (after_restart) { if (after_restart) {
after_restart = false;
if (cvm::proxy->simulation_running()) { if (cvm::proxy->simulation_running()) {
cvm::real const jump2 = dist2(x, x_restart) / (width*width); cvm::real const jump2 = dist2(x, x_restart) / (width*width);
if (jump2 > 0.25) { if (jump2 > 0.25) {
@ -1122,8 +1175,7 @@ int colvar::calc_colvar_properties()
// initialize the restraint center in the first step to the value // initialize the restraint center in the first step to the value
// just calculated from the cvcs // just calculated from the cvcs
// TODO: put it in the restart information if (cvm::step_relative() == 0 && !after_restart) {
if (cvm::step_relative() == 0) {
xr = x; xr = x;
vr.reset(); // (already 0; added for clarity) vr.reset(); // (already 0; added for clarity)
} }
@ -1148,6 +1200,8 @@ int colvar::calc_colvar_properties()
ft_reported = ft; ft_reported = ft;
} }
// At the end of the first update after a restart, we can reset the flag
after_restart = false;
return COLVARS_OK; return COLVARS_OK;
} }
@ -1173,51 +1227,17 @@ cvm::real colvar::update_forces_energy()
f -= fj; 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 // At this point f is the force f from external biases that will be applied to the
// extended variable if there is one // extended variable if there is one
if (is_enabled(f_cv_extended_Lagrangian)) { if (is_enabled(f_cv_extended_Lagrangian)) {
if (cvm::debug()) { 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(); cvm::real dt = cvm::dt();
colvarvalue f_ext(fr.type()); colvarvalue f_ext(fr.type()); // force acting on the extended variable
f_ext.reset(); f_ext.reset();
// the total force is applied to the fictitious mass, while the // 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_ext: total force on extended variable (including harmonic spring)
// f: - initially, external biasing force // f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates // - 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; fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(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); 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 // Now adding the force on the actual colvar (for those biases who
// bypass the extended Lagrangian mass) // bypass the extended Lagrangian mass)
f += fb_actual; f += fb_actual;
@ -1286,8 +1304,10 @@ cvm::real colvar::update_forces_energy()
void colvar::communicate_forces() void colvar::communicate_forces()
{ {
size_t i; size_t i;
if (cvm::debug()) if (cvm::debug()) {
cvm::log("Communicating forces from colvar \""+this->name+"\".\n"); 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)) { if (is_enabled(f_cv_scripted)) {
std::vector<cvm::matrix2d<cvm::real> > func_grads; std::vector<cvm::matrix2d<cvm::real> > func_grads;
@ -1364,7 +1384,7 @@ int colvar::update_cvc_flags()
active_cvc_square_norm = 0.; active_cvc_square_norm = 0.;
for (size_t i = 0; i < cvcs.size(); i++) { 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()) { if (cvcs[i]->is_enabled()) {
n_active_cvcs++; n_active_cvcs++;
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff; active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;

14
lib/colvars/colvar.h
View File

@ -227,11 +227,23 @@ public:
/// Constructor /// Constructor
colvar(std::string const &conf); colvar();
/// Main init function
int init(std::string const &conf);
/// Parse the CVC configuration and allocate their data /// Parse the CVC configuration and allocate their data
int init_components(std::string const &conf); 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: private:
/// Parse the CVC configuration for all components of a certain type /// Parse the CVC configuration for all components of a certain type
template<typename def_class_name> int init_components_type(std::string const &conf, template<typename def_class_name> int init_components_type(std::string const &conf,

14
lib/colvars/colvaratoms.cpp
View File

@ -98,7 +98,7 @@ cvm::atom_group::atom_group()
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); (cvm::proxy)->clear_atom_group(index);
index = -1; 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 // We need to know the fitting options to decide whether the group is scalable
parse_error |= parse_fitting_options(group_conf); 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_com);
enable(f_ag_scalable); 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) int cvm::atom_group::add_index_group(std::string const &index_group_name)
{ {
std::list<std::string>::iterator names_i = cvm::index_group_names.begin(); colvarmodule *cv = cvm::main();
std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) { std::list<std::string>::iterator names_i = cv->index_group_names.begin();
std::list<std::vector<int> >::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) if (*names_i == index_group_name)
break; 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 "+ cvm::error("Error: could not find index group "+
index_group_name+" among those provided by the index file.\n", index_group_name+" among those provided by the index file.\n",
INPUT_ERROR); INPUT_ERROR);

51
lib/colvars/colvarbias.cpp
View File

@ -19,20 +19,6 @@ colvarbias::colvarbias(char const *key)
rank = 1; 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; has_data = false;
b_output_energy = false; b_output_energy = false;
reset(); reset();
@ -48,7 +34,11 @@ int colvarbias::init(std::string const &conf)
colvarparse::init(conf); colvarparse::init(conf);
if (name.size() == 0) { if (name.size() == 0) {
// first initialization
cvm::log("Initializing a new \""+bias_type+"\" instance.\n"); 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)); 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 // lookup the associated colvars
std::vector<std::string> colvar_names; std::vector<std::string> colvar_names;
if (get_keyval(conf, "colvars", 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", cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n",
INPUT_ERROR); INPUT_ERROR);
return 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); cvm::error("Error: no collective variables specified.\n", INPUT_ERROR);
return INPUT_ERROR; return INPUT_ERROR;
} }
@ -89,6 +79,8 @@ int colvarbias::init(std::string const &conf)
cvm::log("Reinitializing bias \""+name+"\".\n"); cvm::log("Reinitializing bias \""+name+"\".\n");
} }
output_prefix = cvm::output_prefix();
get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy); get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
return COLVARS_OK; return COLVARS_OK;
@ -98,7 +90,7 @@ int colvarbias::init(std::string const &conf)
int colvarbias::reset() int colvarbias::reset()
{ {
bias_energy = 0.0; 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(); colvar_forces[i].reset();
} }
return COLVARS_OK; return COLVARS_OK;
@ -132,12 +124,13 @@ int colvarbias::clear()
} }
} }
colvarmodule *cv = cvm::main();
// ...and from the colvars module // ...and from the colvars module
for (std::vector<colvarbias *>::iterator bi = cvm::biases.begin(); for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
bi != cvm::biases.end(); bi != cv->biases.end();
++bi) { ++bi) {
if ( *bi == this) { if ( *bi == this) {
cvm::biases.erase(bi); cv->biases.erase(bi);
break; break;
} }
} }
@ -185,21 +178,29 @@ int colvarbias::add_colvar(std::string const &cv_name)
int colvarbias::update() int colvarbias::update()
{ {
// Note: if anything is added here, it should be added also in the SMP block of calc_biases() if (cvm::debug()) {
// TODO move here debug msg of bias update cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n");
}
has_data = true; has_data = true;
bias_energy = 0.0;
for (size_t ir = 0; ir < num_variables(); ir++) {
colvar_forces[ir].reset();
}
return COLVARS_OK; return COLVARS_OK;
} }
void colvarbias::communicate_forces() 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()) { if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+ 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]);
} }
} }

39
lib/colvars/colvarbias.h
View File

@ -32,17 +32,39 @@ public:
/// Add a new collective variable to this bias /// Add a new collective variable to this bias
int add_colvar(std::string const &cv_name); int add_colvar(std::string const &cv_name);
/// Add a new collective variable to this bias /// How many variables are defined for this bias
size_t number_of_colvars() const inline size_t num_variables() const
{ {
return colvars.size(); return colvars.size();
} }
/// Access the variables vector
inline std::vector<colvar *> *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 /// Retrieve colvar values and calculate their biasing forces
/// Return bias energy /// Return bias energy
virtual int update(); 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<colvarvalue> const &values =
std::vector<colvarvalue>(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 /// Load new configuration - force constant and/or centers only
virtual int change_configuration(std::string const &conf); virtual int change_configuration(std::string const &conf);
@ -51,10 +73,13 @@ public:
virtual cvm::real energy_difference(std::string const &conf); virtual cvm::real energy_difference(std::string const &conf);
/// Give the total number of bins for a given bias. /// Give the total number of bins for a given bias.
// FIXME this is currently 1D only
virtual int bin_num(); virtual int bin_num();
/// Calculate the bin index for a given bias. /// Calculate the bin index for a given bias.
// FIXME this is currently 1D only
virtual int current_bin(); virtual int current_bin();
//// Give the count at a given bin index. //// Give the count at a given bin index.
// FIXME this is currently 1D only
virtual int bin_count(int bin_index); virtual int bin_count(int bin_index);
//// Share information between replicas, whatever it may be. //// Share information between replicas, whatever it may be.
virtual int replica_share(); virtual int replica_share();
@ -62,9 +87,6 @@ public:
/// Perform analysis tasks /// Perform analysis tasks
virtual void analyze() {} virtual void analyze() {}
/// Send forces to the collective variables
virtual void communicate_forces();
/// \brief Constructor /// \brief Constructor
colvarbias(char const *key); colvarbias(char const *key);
@ -135,6 +157,9 @@ public:
return COLVARS_OK; 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 /// If this bias is communicating with other replicas through files, send it to them
virtual int write_state_to_replicas() virtual int write_state_to_replicas()
{ {
@ -162,7 +187,7 @@ protected:
/// through each colvar object /// through each colvar object
std::vector<colvar *> colvars; std::vector<colvar *> colvars;
/// \brief Current forces from this bias to the colvars /// \brief Current forces from this bias to the variables
std::vector<colvarvalue> colvar_forces; std::vector<colvarvalue> colvar_forces;
/// \brief Current energy of this bias (colvar_forces should be obtained by deriving this) /// \brief Current energy of this bias (colvar_forces should be obtained by deriving this)

14
lib/colvars/colvarbias_abf.cpp
View File

@ -30,9 +30,8 @@ int colvarbias_abf::init(std::string const &conf)
{ {
colvarbias::init(conf); colvarbias::init(conf);
provide(f_cvb_scalar_variables);
enable(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 // TODO relax this in case of VMD plugin
if (cvm::temperature() == 0.0) if (cvm::temperature() == 0.0)
@ -221,9 +220,6 @@ colvarbias_abf::~colvarbias_abf()
delete [] system_force; delete [] system_force;
system_force = NULL; 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 // update the output prefix; TODO: move later to setup_output() function
if ( cvm::n_abf_biases == 1 && cvm::n_meta_biases == 0 ) { if (cvm::num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
// This is the only ABF bias // This is the only bias computing PMFs
output_prefix = cvm::output_prefix; output_prefix = cvm::output_prefix();
} else { } else {
output_prefix = cvm::output_prefix + "." + this->name; output_prefix = cvm::output_prefix() + "." + this->name;
} }
if (output_freq && (cvm::step_absolute() % output_freq) == 0) { if (output_freq && (cvm::step_absolute() % output_freq) == 0) {

5
lib/colvars/colvarbias_alb.cpp
View File

@ -40,11 +40,8 @@ int colvarbias_alb::init(std::string const &conf)
{ {
colvarbias::init(conf); colvarbias::init(conf);
provide(f_cvb_scalar_variables);
enable(f_cvb_scalar_variables); enable(f_cvb_scalar_variables);
provide(f_cvb_history_dependent);
size_t i; size_t i;
// get the initial restraint centers // get the initial restraint centers
@ -138,8 +135,6 @@ int colvarbias_alb::init(std::string const &conf)
colvarbias_alb::~colvarbias_alb() colvarbias_alb::~colvarbias_alb()
{ {
if (cvm::n_rest_biases > 0)
cvm::n_rest_biases -= 1;
} }

10
lib/colvars/colvarbias_histogram.cpp
View File

@ -24,10 +24,7 @@ int colvarbias_histogram::init(std::string const &conf)
{ {
colvarbias::init(conf); colvarbias::init(conf);
provide(f_cvb_scalar_variables);
enable(f_cvb_scalar_variables); enable(f_cvb_scalar_variables);
provide(f_cvb_history_dependent);
enable(f_cvb_history_dependent); enable(f_cvb_history_dependent);
size_t i; size_t i;
@ -104,9 +101,6 @@ colvarbias_histogram::~colvarbias_histogram()
delete grid; delete grid;
grid = NULL; 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 // At the first timestep, we need to assign out_name since
// output_prefix is unset during the constructor // output_prefix is unset during the constructor
if (cvm::step_relative() == 0) { 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 + "\""); cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"");
} }
} }
if (out_name_dx.size() == 0) { if (out_name_dx.size() == 0) {
if (cvm::step_relative() == 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 + "\""); cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\"");
} }
} }

473
lib/colvars/colvarbias_meta.cpp
View File

@ -43,7 +43,7 @@ int colvarbias_meta::init(std::string const &conf)
{ {
colvarbias::init(conf); colvarbias::init(conf);
provide(f_cvb_history_dependent); enable(f_cvb_calc_pmf);
get_keyval(conf, "hillWeight", hill_weight, 0.0); get_keyval(conf, "hillWeight", hill_weight, 0.0);
if (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); get_keyval(conf, "hillWidth", hill_width, std::sqrt(2.0 * PI) / 2.0);
cvm::log("Half-widths of the Gaussian hills (sigma's):\n"); cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
cvm::log(colvars[i]->name+std::string(": ")+ cvm::log(variables(i)->name+std::string(": ")+
cvm::to_str(0.5 * colvars[i]->width * hill_width)); 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; comm = single_replica;
} }
// This implies gradients for all colvars // in all cases, the first replica is this bias itself
enable(f_cvb_apply_force); if (replicas.size() == 0) {
replicas.push_back(this);
}
get_keyval(conf, "useGrids", use_grids, true); get_keyval(conf, "useGrids", use_grids, true);
@ -84,14 +86,14 @@ int colvarbias_meta::init(std::string const &conf)
expand_grids = false; expand_grids = false;
size_t i; size_t i;
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
colvars[i]->enable(f_cv_grid); variables(i)->enable(f_cv_grid);
if (colvars[i]->expand_boundaries) { if (variables(i)->expand_boundaries) {
expand_grids = true; expand_grids = true;
cvm::log("Metadynamics bias \""+this->name+"\""+ cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": Will expand grids when the colvar \""+ ": 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); get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) { if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
cvm::log("Option \"saveFreeEnergyFile\" is deprecated, " cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
"please use \"keepFreeEnergyFiles\" instead."); "please use \"keepFreeEnergyFile\" instead.");
} }
get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save); 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); 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 // for well-tempered metadynamics
get_keyval(conf, "wellTempered", well_tempered, false); get_keyval(conf, "wellTempered", well_tempered, false);
get_keyval(conf, "biasTemperature", bias_temperature, -1.0); 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("Well-tempered metadynamics is used.\n");
cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\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 // for ebmeta
target_dist = NULL; target_dist = NULL;
get_keyval(conf, "ebMeta", ebmeta, false); 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); 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; return COLVARS_OK;
} }
@ -234,9 +249,6 @@ colvarbias_meta::~colvarbias_meta()
delete target_dist; delete target_dist;
target_dist = NULL; 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() int colvarbias_meta::update()
{ {
if (cvm::debug()) int error_code = COLVARS_OK;
cvm::log("Updating the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
// 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) { if (use_grids) {
std::vector<int> curr_bin = hills_energy->get_colvars_index(); std::vector<int> 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) { if (expand_grids) {
// first of all, expand the grids, if specified // 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; bool changed_grids = false;
int const min_buffer = int const min_buffer =
(3 * (size_t) std::floor(hill_width)) + 1; (3 * (size_t) std::floor(hill_width)) + 1;
@ -339,9 +373,9 @@ int colvarbias_meta::update()
std::vector<colvarvalue> new_lower_boundaries(hills_energy->lower_boundaries); std::vector<colvarvalue> new_lower_boundaries(hills_energy->lower_boundaries);
std::vector<colvarvalue> new_upper_boundaries(hills_energy->upper_boundaries); std::vector<colvarvalue> 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; continue;
cvm::real &new_lb = new_lower_boundaries[i].real_value; cvm::real &new_lb = new_lower_boundaries[i].real_value;
@ -349,10 +383,10 @@ int colvarbias_meta::update()
int &new_size = new_sizes[i]; int &new_size = new_sizes[i];
bool changed_lb = false, changed_ub = false; 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) { if (curr_bin[i] < min_buffer) {
int const extra_points = (min_buffer - curr_bin[i]); 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; new_size += extra_points;
// changed offset in this direction => the pointer needs to // changed offset in this direction => the pointer needs to
// be changed, too // be changed, too
@ -362,21 +396,21 @@ int colvarbias_meta::update()
cvm::log("Metadynamics bias \""+this->name+"\""+ cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new lower boundary for colvar \""+ ": new lower boundary for colvar \""+
colvars[i]->name+"\", at "+ variables(i)->name+"\", at "+
cvm::to_str(new_lower_boundaries[i])+".\n"); 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) { if (curr_bin[i] > new_size - min_buffer - 1) {
int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer); 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; new_size += extra_points;
changed_ub = true; changed_ub = true;
cvm::log("Metadynamics bias \""+this->name+"\""+ cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new upper boundary for colvar \""+ ": new upper boundary for colvar \""+
colvars[i]->name+"\", at "+ variables(i)->name+"\", at "+
cvm::to_str(new_upper_boundaries[i])+".\n"); 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->lower_boundaries = new_lower_boundaries;
new_hills_energy_gradients->upper_boundaries = new_upper_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->map_grid(*hills_energy);
new_hills_energy_gradients->map_grid(*hills_energy_gradients); 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 // 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+"\""+ cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n"); ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n");
}
cvm::real hills_scale=1.0; cvm::real hills_scale=1.0;
if (ebmeta) { if (ebmeta) {
hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index()); hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
if(cvm::step_absolute() <= long(ebmeta_equil_steps)) { if(cvm::step_absolute() <= long(ebmeta_equil_steps)) {
cvm::real const hills_lambda = cvm::real const hills_lambda =
(cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) / (cvm::real(long(ebmeta_equil_steps) - cvm::step_absolute())) /
(cvm::real(ebmeta_equil_steps)); (cvm::real(ebmeta_equil_steps));
hills_scale = hills_lambda + (1-hills_lambda)*hills_scale; hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
} }
} }
if (well_tempered) { if (well_tempered) {
@ -471,160 +512,165 @@ int colvarbias_meta::update()
} }
} }
// sync with the other replicas (if needed) return COLVARS_OK;
if (comm != single_replica) { }
// 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 return COLVARS_OK;
bias_energy = 0.0; }
for (size_t ir = 0; ir < colvars.size(); ir++) {
colvar_forces[ir].reset();
int colvarbias_meta::calc_energy(std::vector<colvarvalue> 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++) { std::vector<int> const curr_bin = values.size() ?
replicas[ir]->bias_energy = 0.0; hills_energy->get_colvars_index(values) :
for (size_t ic = 0; ic < colvars.size(); ic++) { hills_energy->get_colvars_index();
replicas[ir]->colvar_forces[ic].reset();
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");
} }
} }
} else {
if (use_grids) { // off the grid: compute analytically only the hills at the grid's edges
for (ir = 0; ir < replicas.size(); ir++) {
// get the forces from the grid calc_hills(replicas[ir]->hills_off_grid.begin(),
replicas[ir]->hills_off_grid.end(),
if ((cvm::step_absolute() % grids_freq) == 0) { bias_energy,
// map the most recent gaussians to the grids values);
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<int> curr_bin = hills_energy->get_colvars_index();
if (cvm::debug())
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": current coordinates on the grid: "+
cvm::to_str(curr_bin)+".\n");
if (hills_energy->index_ok(curr_bin)) {
// within the grid: add the energy and the forces from there
bias_energy += hills_energy->value(curr_bin);
for (size_t ic = 0; ic < colvars.size(); ic++) {
cvm::real const *f = &(hills_energy_gradients->value(curr_bin));
colvar_forces[ic].real_value += -1.0 * f[ic];
// the gradients are stored, not the forces
}
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
bias_energy += replicas[ir]->hills_energy->value(curr_bin);
cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
for (size_t ic = 0; ic < colvars.size(); ic++) {
colvar_forces[ic].real_value += -1.0 * f[ic];
}
}
} else {
// off the grid: compute analytically only the hills at the grid's edges
calc_hills(hills_off_grid.begin(), hills_off_grid.end(), bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force(ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
}
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
calc_hills(replicas[ir]->hills_off_grid.begin(),
replicas[ir]->hills_off_grid.end(),
bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force(ic,
replicas[ir]->hills_off_grid.begin(),
replicas[ir]->hills_off_grid.end(),
colvar_forces);
}
}
} }
} }
// now include the hills that have not been binned yet (starting // now include the hills that have not been binned yet (starting
// from new_hills_begin) // from new_hills_begin)
calc_hills(new_hills_begin, hills.end(), bias_energy); for (ir = 0; ir < replicas.size(); ir++) {
for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills(replicas[ir]->new_hills_begin,
calc_hills_force(ic, new_hills_begin, hills.end(), colvar_forces); replicas[ir]->hills.end(),
} bias_energy);
if (cvm::debug()) {
if (cvm::debug()) cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n");
cvm::log("Hills energy = "+cvm::to_str(bias_energy)+
", hills forces = "+cvm::to_str(colvar_forces)+".\n");
if (cvm::debug())
cvm::log("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding the forces from the other replicas.\n");
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
calc_hills(replicas[ir]->new_hills_begin,
replicas[ir]->hills.end(),
bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force(ic,
replicas[ir]->new_hills_begin,
replicas[ir]->hills.end(),
colvar_forces);
}
if (cvm::debug())
cvm::log("Hills energy = "+cvm::to_str(bias_energy)+
", hills forces = "+cvm::to_str(colvar_forces)+".\n");
} }
}
return COLVARS_OK; return COLVARS_OK;
} }
int colvarbias_meta::calc_forces(std::vector<colvarvalue> 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<int> 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, void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last, colvarbias_meta::hill_iter h_last,
cvm::real &energy, cvm::real &energy,
std::vector<colvarvalue> const &colvar_values) std::vector<colvarvalue> const &colvar_values)
{ {
std::vector<colvarvalue> curr_values(colvars.size()); size_t i = 0;
for (size_t i = 0; i < colvars.size(); i++) { std::vector<colvarvalue> curr_values(num_variables());
curr_values[i].type(colvars[i]->value()); for (i = 0; i < num_variables(); i++) {
curr_values[i].type(variables(i)->value());
} }
if (colvar_values.size()) { 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]; curr_values[i] = colvar_values[i];
} }
} else { } else {
for (size_t i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
curr_values[i] = colvars[i]->value(); 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 // compute the gaussian exponent
cvm::real cv_sqdev = 0.0; 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 &x = curr_values[i];
colvarvalue const &center = h->centers[i]; colvarvalue const &center = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[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 // compute the gaussian
@ -658,14 +704,14 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
std::vector<colvarvalue> const &values) std::vector<colvarvalue> const &values)
{ {
// Retrieve the value of the colvar // 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 // do the type check only once (all colvarvalues in the hills series
// were already saved with their types matching those in the // were already saved with their types matching those in the
// colvars) // colvars)
hill_iter h; hill_iter h;
switch (colvars[i]->value().type()) { switch (variables(i)->value().type()) {
case colvarvalue::type_scalar: case colvarvalue::type_scalar:
for (h = h_first; h != h_last; h++) { 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]; cvm::real const half_width = 0.5 * h->widths[i];
forces[i].real_value += forces[i].real_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) * ( 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; 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]; cvm::real const half_width = 0.5 * h->widths[i];
forces[i].rvector_value += forces[i].rvector_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) * ( 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; 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]; cvm::real const half_width = 0.5 * h->widths[i];
forces[i].quaternion_value += forces[i].quaternion_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) * ( 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; 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]; cvm::real const half_width = 0.5 * h->widths[i];
forces[i].vector1d_value += forces[i].vector1d_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) * ( 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; 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 // TODO: improve it by looping over a small subgrid instead of the whole grid
std::vector<colvarvalue> colvar_values(colvars.size()); std::vector<colvarvalue> colvar_values(num_variables());
std::vector<cvm::real> colvar_forces_scalar(colvars.size()); std::vector<cvm::real> colvar_forces_scalar(num_variables());
std::vector<int> he_ix = he->new_index(); std::vector<int> he_ix = he->new_index();
std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0); std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
cvm::real hills_energy_here = 0.0; cvm::real hills_energy_here = 0.0;
std::vector<colvarvalue> hills_forces_here(colvars.size(), 0.0); std::vector<colvarvalue> hills_forces_here(num_variables(), 0.0);
size_t count = 0; size_t count = 0;
size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1))); 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)); (he->index_ok(he_ix)) && (hg->index_ok(hg_ix));
count++) { count++) {
size_t i; 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); 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); calc_hills(h_first, h_last, hills_energy_here, colvar_values);
he->acc_value(he_ix, hills_energy_here); 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(); hills_forces_here[i].reset();
calc_hills_force(i, h_first, h_last, hills_forces_here, colvar_values); calc_hills_force(i, h_first, h_last, hills_forces_here, colvar_values);
colvar_forces_scalar[i] = hills_forces_here[i].real_value; 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 ( ; (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); 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() void colvarbias_meta::update_replicas_registry()
{ {
if (cvm::debug()) if (cvm::debug())
@ -975,7 +1036,6 @@ void colvarbias_meta::update_replicas_registry()
} }
} }
if (cvm::debug()) if (cvm::debug())
cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+ cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
cvm::to_str(replicas.size())+" elements.\n"); cvm::to_str(replicas.size())+" elements.\n");
@ -984,7 +1044,8 @@ void colvarbias_meta::update_replicas_registry()
void colvarbias_meta::read_replica_files() 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 (! (replicas[ir])->replica_state_file_in_sync) {
// if a new state file is being read, the hills file is also new // 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; cvm::real h_weight;
get_keyval(data, "weight", h_weight, hill_weight, parse_silent); get_keyval(data, "weight", h_weight, hill_weight, parse_silent);
std::vector<colvarvalue> h_centers(colvars.size()); std::vector<colvarvalue> h_centers(num_variables());
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
h_centers[i].type(colvars[i]->value()); h_centers[i].type(variables(i)->value());
} }
{ {
// it is safer to read colvarvalue objects one at a time; // 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; std::string centers_input;
key_lookup(data, "centers", centers_input); key_lookup(data, "centers", centers_input);
std::istringstream centers_is(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]; centers_is >> h_centers[i];
} }
} }
std::vector<cvm::real> h_widths(colvars.size()); std::vector<cvm::real> h_widths(num_variables());
get_keyval(data, "widths", h_widths, get_keyval(data, "widths", h_widths,
std::vector<cvm::real> (colvars.size(), (std::sqrt(2.0 * PI) / 2.0)), std::vector<cvm::real>(num_variables(), (std::sqrt(2.0 * PI) / 2.0)),
parse_silent); parse_silent);
std::string h_replica = ""; std::string h_replica = "";
@ -1406,6 +1467,13 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
int colvarbias_meta::setup_output() 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) { if (comm == multiple_replicas) {
@ -1421,10 +1489,10 @@ int colvarbias_meta::setup_output()
// those by another replica // those by another replica
replica_hills_file = replica_hills_file =
(std::string(pwd)+std::string(PATHSEP)+ (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 = replica_state_file =
(std::string(pwd)+std::string(PATHSEP)+ (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; delete[] pwd;
// now register this replica // now register this replica
@ -1492,13 +1560,14 @@ int colvarbias_meta::setup_output()
} }
if (b_hills_traj) { 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+ ".colvars."+this->name+
( (comm != single_replica) ? ( (comm != single_replica) ?
("."+replica_id) : ("."+replica_id) :
("") )+ ("") )+
".hills.traj"); ".hills.traj");
if (!hills_traj_os.is_open()) { if (!hills_traj_os.is_open()) {
cvm::backup_file(traj_file_name.c_str());
hills_traj_os.open(traj_file_name.c_str()); hills_traj_os.open(traj_file_name.c_str());
} }
if (!hills_traj_os.is_open()) 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); colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy);
pmf->setup(); 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)) { if ((comm == single_replica) || (dump_replica_fes)) {
// output the PMF from this instance or replica // output the PMF from this instance or replica
pmf->reset(); pmf->reset();
@ -1607,7 +1666,7 @@ void colvarbias_meta::write_pmf()
pmf->multiply_constant(well_temper_scale); 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" : "") + ((comm != single_replica) ? ".partial" : "") +
(dump_fes_save ? (dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") + "."+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; cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
pmf->multiply_constant(well_temper_scale); 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 ? (dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") + "."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf"); ".pmf");

26
lib/colvars/colvarbias_meta.h
View File

@ -36,8 +36,20 @@ public:
colvarbias_meta(char const *key); colvarbias_meta(char const *key);
virtual int init(std::string const &conf); 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 ~colvarbias_meta();
virtual int update(); 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<colvarvalue> const &values =
std::vector<colvarvalue>(0));
virtual int calc_forces(std::vector<colvarvalue> const &values =
std::vector<colvarvalue>(0));
virtual std::string const get_state_params() const; virtual std::string const get_state_params() const;
virtual int set_state_params(std::string const &state_conf); virtual int set_state_params(std::string const &state_conf);
@ -102,18 +114,18 @@ protected:
/// \brief Calculate the values of the hills, incrementing /// \brief Calculate the values of the hills, incrementing
/// bias_energy /// bias_energy
virtual void calc_hills(hill_iter h_first, virtual void calc_hills(hill_iter h_first,
hill_iter h_last, hill_iter h_last,
cvm::real &energy, cvm::real &energy,
std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0)); std::vector<colvarvalue> const &values = std::vector<colvarvalue>(0));
/// \brief Calculate the forces acting on the i-th colvar, /// \brief Calculate the forces acting on the i-th colvar,
/// incrementing colvar_forces[i]; must be called after calc_hills /// incrementing colvar_forces[i]; must be called after calc_hills
/// each time the values of the colvars are changed /// each time the values of the colvars are changed
virtual void calc_hills_force(size_t const &i, virtual void calc_hills_force(size_t const &i,
hill_iter h_first, hill_iter h_first,
hill_iter h_last, hill_iter h_last,
std::vector<colvarvalue> &forces, std::vector<colvarvalue> &forces,
std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0)); std::vector<colvarvalue> const &values = std::vector<colvarvalue>(0));
/// Height of new hills /// Height of new hills

214
lib/colvars/colvarbias_restraint.cpp
View File

@ -33,17 +33,15 @@ int colvarbias_restraint::init(std::string const &conf)
int colvarbias_restraint::update() int colvarbias_restraint::update()
{ {
bias_energy = 0.0; // Update base class (bias_energy and colvar_forces are zeroed there)
colvarbias::update();
if (cvm::debug())
cvm::log("Updating the restraint bias \""+this->name+"\".\n");
// Force and energy calculation // Force and energy calculation
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
colvar_forces[i].type(colvars[i]->value()); bias_energy += restraint_potential(i);
colvar_forces[i].type(variables(i)->value());
colvar_forces[i].is_derivative(); colvar_forces[i].is_derivative();
colvar_forces[i] = restraint_force(i); colvar_forces[i] = restraint_force(i);
bias_energy += restraint_potential(i);
} }
if (cvm::debug()) if (cvm::debug())
@ -59,8 +57,6 @@ int colvarbias_restraint::update()
colvarbias_restraint::~colvarbias_restraint() 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); bool null_centers = (colvar_centers.size() == 0);
if (null_centers) { if (null_centers) {
// try to initialize the restraint centers for the first time // try to initialize the restraint centers for the first time
colvar_centers.resize(colvars.size()); colvar_centers.resize(num_variables());
colvar_centers_raw.resize(colvars.size()); colvar_centers_raw.resize(num_variables());
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
colvar_centers[i].type(colvars[i]->value()); colvar_centers[i].type(variables(i)->value());
colvar_centers[i].reset(); 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(); colvar_centers_raw[i].reset();
} }
} }
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { 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()) { if (cvm::debug()) {
cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n"); 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; return INPUT_ERROR;
} }
if (colvar_centers.size() != colvars.size()) { if (colvar_centers.size() != num_variables()) {
cvm::error("Error: number of centers does not match " cvm::error("Error: number of centers does not match "
"that of collective variables.\n", INPUT_ERROR); "that of collective variables.\n", INPUT_ERROR);
return 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) int colvarbias_restraint_centers::change_configuration(std::string const &conf)
{ {
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
colvar_centers[i].type(colvars[i]->value()); colvar_centers[i].type(variables(i)->value());
colvar_centers[i].apply_constraints(); 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]; colvar_centers_raw[i] = colvar_centers[i];
} }
} }
@ -269,7 +265,7 @@ int colvarbias_restraint_centers_moving::init(std::string const &conf)
size_t i; size_t i;
if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { 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 " cvm::error("Error: number of target centers does not match "
"that of collective variables.\n"); "that of collective variables.\n");
} }
@ -308,9 +304,9 @@ int colvarbias_restraint_centers_moving::update()
// at each simulation step (or stage, if applicable) // at each simulation step (or stage, if applicable)
// (take current stage into account: it can be non-zero // (take current stage into account: it can be non-zero
// if we are restarting a staged calculation) // if we are restarting a staged calculation)
centers_incr.resize(colvars.size()); centers_incr.resize(num_variables());
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
centers_incr[i].type(colvars[i]->value()); centers_incr[i].type(variables(i)->value());
centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) /
cvm::real( target_nstages ? (target_nstages - stage) : cvm::real( target_nstages ? (target_nstages - stage) :
(target_nsteps - cvm::step_absolute())); (target_nsteps - cvm::step_absolute()));
@ -326,10 +322,10 @@ int colvarbias_restraint_centers_moving::update()
&& (cvm::step_absolute() % target_nsteps) == 0 && (cvm::step_absolute() % target_nsteps) == 0
&& stage < target_nstages) { && 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_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[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(); colvar_centers[i].apply_constraints();
} }
stage++; stage++;
@ -341,10 +337,10 @@ int colvarbias_restraint_centers_moving::update()
} else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) { } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) {
// move the restraint centers in the direction of the targets // move the restraint centers in the direction of the targets
// (slow growth) // (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_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[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(); colvar_centers[i].apply_constraints();
} }
} }
@ -363,7 +359,7 @@ int colvarbias_restraint_centers_moving::update_acc_work()
{ {
if (b_output_acc_work) { if (b_output_acc_work) {
if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { 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 // project forces on the calculated increments at this step
acc_work += colvar_forces[i] * centers_incr[i]; 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) { if (b_chg_centers) {
size_t i; size_t i;
os << "centers "; os << "centers ";
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
os << " " os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers[i]; << colvar_centers[i];
} }
os << "\n"; os << "\n";
os << "centers_raw "; os << "centers_raw ";
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
os << " " os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers_raw[i]; << 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) std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os)
{ {
if (b_output_centers) { if (b_output_centers) {
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
size_t const this_cv_width = (colvars[i]->value()).output_width(cvm::cv_width); size_t const this_cv_width = (variables(i)->value()).output_width(cvm::cv_width);
os << " x0_" 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) std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
{ {
if (b_output_centers) { if (b_output_centers) {
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
os << " " os << " "
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
<< colvar_centers[i]; << colvar_centers[i];
@ -539,9 +535,9 @@ int colvarbias_restraint_k_moving::update()
} }
force_k = starting_force_k + (target_force_k - starting_force_k) force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow(lambda, force_k_exp); * std::pow(lambda, force_k_exp);
cvm::log("Restraint " + this->name + ", stage " + cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); + " : lambda = " + cvm::to_str(lambda)
cvm::log("Setting force constant to " + cvm::to_str(force_k)); + ", k = " + cvm::to_str(force_k));
} }
// TI calculation: estimate free energy derivative // TI calculation: estimate free energy derivative
@ -557,7 +553,7 @@ int colvarbias_restraint_k_moving::update()
// Square distance normalized by square colvar width // Square distance normalized by square colvar width
cvm::real dist_sq = 0.0; 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); dist_sq += d_restraint_potential_dk(i);
} }
@ -569,7 +565,8 @@ int colvarbias_restraint_k_moving::update()
if (cvm::step_absolute() % target_nsteps == 0 && if (cvm::step_absolute() % target_nsteps == 0 &&
cvm::step_absolute() > 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))); + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps)));
// ...and move on to the next one // ...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) force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow(lambda, force_k_exp); * std::pow(lambda, force_k_exp);
cvm::log("Restraint " + this->name + ", stage " + cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); + " : lambda = " + cvm::to_str(lambda)
cvm::log("Setting force constant to " + cvm::to_str(force_k)); + ", 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_centers_moving::init(conf);
colvarbias_restraint_k_moving::init(conf); colvarbias_restraint_k_moving::init(conf);
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
if (colvars[i]->width != 1.0) if (variables(i)->width != 1.0)
cvm::log("The force constant for colvar \""+colvars[i]->name+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+ "\" 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"); " 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 cvm::real colvarbias_restraint_harmonic::restraint_potential(size_t i) const
{ {
return 0.5 * force_k / (colvars[i]->width * colvars[i]->width) * return 0.5 * force_k / (variables(i)->width * variables(i)->width) *
colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]); variables(i)->dist2(variables(i)->value(), colvar_centers[i]);
} }
colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const
{ {
return -0.5 * force_k / (colvars[i]->width * colvars[i]->width) * return -0.5 * force_k / (variables(i)->width * variables(i)->width) *
colvars[i]->dist2_lgrad(colvars[i]->value(), colvar_centers[i]); variables(i)->dist2_lgrad(variables(i)->value(), colvar_centers[i]);
} }
cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const
{ {
return 0.5 / (colvars[i]->width * colvars[i]->width) * return 0.5 / (variables(i)->width * variables(i)->width) *
colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]); 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_moving(key),
colvarbias_restraint_k_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_moving::init(conf);
colvarbias_restraint_k_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); enable(f_cvb_scalar_variables);
size_t i; size_t i;
@ -857,9 +860,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
bool b_null_lower_walls = false; bool b_null_lower_walls = false;
if (lower_walls.size() == 0) { if (lower_walls.size() == 0) {
b_null_lower_walls = true; b_null_lower_walls = true;
lower_walls.resize(number_of_colvars()); lower_walls.resize(num_variables());
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
lower_walls[i].type(colvars[i]->value()); lower_walls[i].type(variables(i)->value());
lower_walls[i].reset(); lower_walls[i].reset();
} }
} }
@ -872,9 +875,9 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
bool b_null_upper_walls = false; bool b_null_upper_walls = false;
if (upper_walls.size() == 0) { if (upper_walls.size() == 0) {
b_null_upper_walls = true; b_null_upper_walls = true;
upper_walls.resize(number_of_colvars()); upper_walls.resize(num_variables());
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
upper_walls[i].type(colvars[i]->value()); upper_walls[i].type(variables(i)->value());
upper_walls[i].reset(); 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)) { if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) {
for (i = 0; i < colvars.size(); i++) { for (i = 0; i < num_variables(); i++) {
if (colvars[i]->is_enabled(f_cv_periodic)) { if (variables(i)->is_enabled(f_cv_periodic)) {
cvm::error("Error: at least one variable is 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; return INPUT_ERROR;
} }
} }
} }
if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) { 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]) { if (lower_walls[i] >= upper_walls[i]) {
cvm::error("Error: one upper wall, "+ cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+ cvm::to_str(upper_walls[i])+
@ -909,13 +912,24 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
INPUT_ERROR); 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++) { for (i = 0; i < num_variables(); i++) {
if (colvars[i]->width != 1.0) if (variables(i)->width != 1.0)
cvm::log("The force constant for colvar \""+colvars[i]->name+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+ "\" 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"); " according to the specified width.\n");
} }
@ -935,20 +949,20 @@ int colvarbias_restraint_harmonic_walls::update()
void colvarbias_restraint_harmonic_walls::communicate_forces() 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()) { if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+ 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 cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
{ {
colvar *cv = colvars[i]; colvar *cv = variables(i);
colvarvalue const &cvv = colvars[i]->actual_value(); colvarvalue const &cvv = variables(i)->actual_value();
// For a periodic colvar, both walls may be applicable at the same time // For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one // 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]); cvm::real const upper_wall_dist2 = cv->dist2(cvv, upper_walls[i]);
if (lower_wall_dist2 < upper_wall_dist2) { if (lower_wall_dist2 < upper_wall_dist2) {
cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); 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 { } else {
cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); 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; return 0.0;
} }
if (lower_walls.size() > 0) { if (lower_walls.size() > 0) {
cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); 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) { if (upper_walls.size() > 0) {
cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); 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; 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 colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) const
{ {
cvm::real const dist = colvar_distance(i); 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; 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 colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const
{ {
cvm::real const dist = colvar_distance(i); 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;
dist; 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 colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i) const
{ {
cvm::real const dist = colvar_distance(i); 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; dist * dist;
} }
@ -1054,16 +1070,16 @@ int colvarbias_restraint_linear::init(std::string const &conf)
colvarbias_restraint_centers_moving::init(conf); colvarbias_restraint_centers_moving::init(conf);
colvarbias_restraint_k_moving::init(conf); colvarbias_restraint_k_moving::init(conf);
for (size_t i = 0; i < colvars.size(); i++) { for (size_t i = 0; i < num_variables(); i++) {
if (colvars[i]->is_enabled(f_cv_periodic)) { if (variables(i)->is_enabled(f_cv_periodic)) {
cvm::error("Error: linear biases cannot be applied to periodic variables.\n", cvm::error("Error: linear biases cannot be applied to periodic variables.\n",
INPUT_ERROR); INPUT_ERROR);
return INPUT_ERROR; return INPUT_ERROR;
} }
if (colvars[i]->width != 1.0) if (variables(i)->width != 1.0)
cvm::log("The force constant for colvar \""+colvars[i]->name+ cvm::log("The force constant for colvar \""+variables(i)->name+
"\" will be rescaled to "+ "\" 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"); " 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 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 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 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 vector_size = 0;
size_t icv; size_t icv;
for (icv = 0; icv < colvars.size(); icv++) { for (icv = 0; icv < num_variables(); icv++) {
vector_size += colvars[icv]->value().size(); vector_size += variables(icv)->value().size();
} }
cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size); cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
// calculate the histogram // calculate the histogram
p.reset(); p.reset();
for (icv = 0; icv < colvars.size(); icv++) { for (icv = 0; icv < num_variables(); icv++) {
colvarvalue const &cv = colvars[icv]->value(); colvarvalue const &cv = variables(icv)->value();
if (cv.type() == colvarvalue::type_scalar) { if (cv.type() == colvarvalue::type_scalar) {
cvm::real const cv_value = cv.real_value; cvm::real const cv_value = cv.real_value;
size_t igrid; size_t igrid;
@ -1309,7 +1325,9 @@ int colvarbias_restraint_histogram::update()
} }
} }
} else { } 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; bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
// calculate the forces // calculate the forces
for (icv = 0; icv < colvars.size(); icv++) { for (icv = 0; icv < num_variables(); icv++) {
colvarvalue const &cv = colvars[icv]->value(); colvarvalue const &cv = variables(icv)->value();
colvarvalue &cv_force = colvar_forces[icv]; colvarvalue &cv_force = colvar_forces[icv];
cv_force.type(cv); cv_force.type(cv);
cv_force.reset(); cv_force.reset();
@ -1363,10 +1381,10 @@ int colvarbias_restraint_histogram::update()
std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
{ {
if (b_write_histogram) { 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()); std::ofstream os(file_name.c_str());
os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width) os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
<< " " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3) << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3)
<< ")\n"; << ")\n";
size_t igrid; size_t igrid;
for (igrid = 0; igrid < p.size(); igrid++) { for (igrid = 0; igrid < p.size(); igrid++) {

6
lib/colvars/colvarbias_restraint.h
View File

@ -260,6 +260,12 @@ protected:
/// \brief Location of the upper walls /// \brief Location of the upper walls
std::vector<colvarvalue> upper_walls; std::vector<colvarvalue> 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 colvar_distance(size_t i) const;
virtual cvm::real restraint_potential(size_t i) const; virtual cvm::real restraint_potential(size_t i) const;
virtual colvarvalue const restraint_force(size_t i) const; virtual colvarvalue const restraint_force(size_t i) const;

20
lib/colvars/colvarcomp.cpp
View File

@ -48,8 +48,6 @@ colvar::cvc::cvc(std::string const &conf)
get_keyval(conf, "period", period, 0.0); get_keyval(conf, "period", period, 0.0);
get_keyval(conf, "wrapAround", wrap_center, 0.0); get_keyval(conf, "wrapAround", wrap_center, 0.0);
// All cvcs implement this
provide(f_cvc_debug_gradient);
get_keyval_feature((colvarparse *)this, conf, "debugGradients", get_keyval_feature((colvarparse *)this, conf, "debugGradients",
f_cvc_debug_gradient, false, parse_silent); 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) 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", if (get_keyval_feature(this, conf, "oneSiteSystemForce",
f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) { f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: " 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))) { f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
cvm::log("Computing total force on group 1 only"); 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<cvm::atom_group *>::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; return COLVARS_OK;
} }
@ -87,8 +100,7 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
group->key = group_key; group->key = group_key;
if (b_try_scalable) { if (b_try_scalable) {
// TODO rewrite this logic in terms of dependencies if (is_available(f_cvc_scalable_com) && is_enabled(f_cvc_com_based)) {
if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
enable(f_cvc_scalable_com); enable(f_cvc_scalable_com);
enable(f_cvc_scalable); enable(f_cvc_scalable);
// The CVC makes the feature available; // The CVC makes the feature available;

9
lib/colvars/colvarcomp.h
View File

@ -343,6 +343,15 @@ public:
virtual void calc_value(); virtual void calc_value();
virtual void calc_gradients(); virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force); 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;
}; };

8
lib/colvars/colvarcomp_angles.cpp
View File

@ -21,7 +21,7 @@ colvar::angle::angle(std::string const &conf)
function_type = "angle"; function_type = "angle";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
group1 = parse_group(conf, "group1"); group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2"); group2 = parse_group(conf, "group2");
@ -40,7 +40,7 @@ colvar::angle::angle(cvm::atom const &a1,
function_type = "angle"; function_type = "angle";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1)); group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2)); group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@ -259,7 +259,7 @@ colvar::dihedral::dihedral(std::string const &conf)
b_periodic = true; b_periodic = true;
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
group1 = parse_group(conf, "group1"); group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2"); group2 = parse_group(conf, "group2");
@ -285,7 +285,7 @@ colvar::dihedral::dihedral(cvm::atom const &a1,
b_periodic = true; b_periodic = true;
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
b_1site_force = false; b_1site_force = false;

41
lib/colvars/colvarcomp_distances.cpp
View File

@ -23,7 +23,7 @@ colvar::distance::distance(std::string const &conf)
function_type = "distance"; function_type = "distance";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
group1 = parse_group(conf, "group1"); group1 = parse_group(conf, "group1");
group2 = parse_group(conf, "group2"); group2 = parse_group(conf, "group2");
@ -44,7 +44,7 @@ colvar::distance::distance()
function_type = "distance"; function_type = "distance";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
b_no_PBC = false; b_no_PBC = false;
x.type(colvarvalue::type_scalar); x.type(colvarvalue::type_scalar);
} }
@ -106,7 +106,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
: distance(conf) : distance(conf)
{ {
function_type = "distance_vec"; function_type = "distance_vec";
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_3vector); x.type(colvarvalue::type_3vector);
} }
@ -115,7 +115,7 @@ colvar::distance_vec::distance_vec()
: distance() : distance()
{ {
function_type = "distance_vec"; function_type = "distance_vec";
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_3vector); x.type(colvarvalue::type_3vector);
} }
@ -176,7 +176,7 @@ colvar::distance_z::distance_z(std::string const &conf)
function_type = "distance_z"; function_type = "distance_z";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar); x.type(colvarvalue::type_scalar);
// TODO detect PBC from MD engine (in simple cases) // TODO detect PBC from MD engine (in simple cases)
@ -228,7 +228,7 @@ colvar::distance_z::distance_z()
function_type = "distance_z"; function_type = "distance_z";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar); x.type(colvarvalue::type_scalar);
} }
@ -372,7 +372,7 @@ colvar::distance_xy::distance_xy(std::string const &conf)
function_type = "distance_xy"; function_type = "distance_xy";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar); x.type(colvarvalue::type_scalar);
} }
@ -383,7 +383,7 @@ colvar::distance_xy::distance_xy()
function_type = "distance_xy"; function_type = "distance_xy";
provide(f_cvc_inv_gradient); provide(f_cvc_inv_gradient);
provide(f_cvc_Jacobian); provide(f_cvc_Jacobian);
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_scalar); x.type(colvarvalue::type_scalar);
} }
@ -479,7 +479,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
: distance(conf) : distance(conf)
{ {
function_type = "distance_dir"; function_type = "distance_dir";
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_unit3vector); x.type(colvarvalue::type_unit3vector);
} }
@ -488,7 +488,7 @@ colvar::distance_dir::distance_dir()
: distance() : distance()
{ {
function_type = "distance_dir"; function_type = "distance_dir";
provide(f_cvc_com_based); enable(f_cvc_com_based);
x.type(colvarvalue::type_unit3vector); 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) colvar::distance_inv::distance_inv(std::string const &conf)
: distance(conf) : distance(conf)

272
lib/colvars/colvardeps.cpp
View File

@ -14,9 +14,6 @@
colvardeps::~colvardeps() { colvardeps::~colvardeps() {
size_t i; size_t i;
for (i=0; i<feature_states.size(); i++) {
if (feature_states[i] != NULL) delete feature_states[i];
}
// Do not delete features if it's static // Do not delete features if it's static
// for (i=0; i<features.size(); i++) { // for (i=0; i<features.size(); i++) {
// if (features[i] != NULL) delete features[i]; // if (features[i] != NULL) delete features[i];
@ -34,16 +31,34 @@ colvardeps::~colvardeps() {
} }
void colvardeps::provide(int feature_id) { void colvardeps::provide(int feature_id, bool truefalse) {
feature_states[feature_id]->available = true; 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, bool colvardeps::get_keyval_feature(colvarparse *cvp,
std::string const &conf, char const *key, std::string const &conf, char const *key,
int feature_id, bool const &def_value, int feature_id, bool const &def_value,
colvarparse::Parse_Mode const parse_mode) 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 value;
bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode); bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode);
if (value) enable(feature_id); if (value) enable(feature_id);
@ -52,19 +67,19 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
int colvardeps::enable(int feature_id, int colvardeps::enable(int feature_id,
bool dry_run /* default: false */, bool dry_run /* default: false */,
// dry_run: fail silently, do not enable if available // dry_run: fail silently, do not enable if available
// flag is passed recursively to deps of this feature // flag is passed recursively to deps of this feature
bool toplevel /* default: true */) bool toplevel /* default: true */)
// toplevel: false if this is called as part of a chain of dependency resolution // 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 // this is used to diagnose failed dependencies by displaying the full stack
// only the toplevel dependency will throw a fatal error // only the toplevel dependency will throw a fatal error
{ {
int res; int res;
size_t i, j; size_t i, j;
bool ok; bool ok;
feature *f = features()[feature_id]; feature *f = features()[feature_id];
feature_state *fs = feature_states[feature_id]; feature_state *fs = &feature_states[feature_id];
if (cvm::debug()) { if (cvm::debug()) {
cvm::log("DEPS: " + description + cvm::log("DEPS: " + description +
@ -88,6 +103,14 @@ int colvardeps::enable(int feature_id,
return COLVARS_ERROR; 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 // 1) enforce exclusions
for (i=0; i<f->requires_exclude.size(); i++) { for (i=0; i<f->requires_exclude.size(); i++) {
feature *g = features()[f->requires_exclude[i]]; feature *g = features()[f->requires_exclude[i]];
@ -168,9 +191,9 @@ int colvardeps::enable(int feature_id,
if (res != COLVARS_OK) { if (res != COLVARS_OK) {
if (!dry_run) { if (!dry_run) {
cvm::log("...required by \"" + f->description + "\" in " + description); cvm::log("...required by \"" + f->description + "\" in " + description);
} if (toplevel) {
if (toplevel) { cvm::error("Error: Failed dependency in " + description + ".");
cvm::error("Error: Failed dependency in " + description + "."); }
} }
return res; return res;
} }
@ -194,9 +217,12 @@ int colvardeps::enable(int feature_id,
// // we need refs to parents to walk up the deps tree! // // we need refs to parents to walk up the deps tree!
// // or refresh // // 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 // Shorthand macros for describing dependencies
#define f_description(f, d) features()[f]->description = d
#define f_req_self(f, g) features()[f]->requires_self.push_back(g) #define f_req_self(f, g) features()[f]->requires_self.push_back(g)
// This macro ensures that exclusions are symmetric // This macro ensures that exclusions are symmetric
#define f_req_exclude(f, g) features()[f]->requires_exclude.push_back(g); \ #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++) { for (i = 0; i < f_cvb_ntot; i++) {
features().push_back(new feature); 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 // Initialize feature_states for each instance
feature_states.reserve(f_cvb_ntot); feature_states.reserve(f_cvb_ntot);
for (i = 0; i < f_cvb_ntot; i++) { 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 // Most features are available, so we set them so
// and list exceptions below // 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); 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); f_req_children(f_cv_active, f_cvc_active);
// Colvars must be either a linear combination, or scalar (and polynomial) or scripted // 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_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_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_gradient);
f_req_self(f_cv_collect_gradient, f_cv_scalar); 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 // 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); f_req_alt2(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
// Deps for explicit total force calculation // 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_scalar);
f_req_self(f_cv_total_force_calc, f_cv_linear); 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_children(f_cv_total_force_calc, f_cvc_inv_gradient);
f_req_self(f_cv_total_force_calc, f_cv_Jacobian); 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_scalar);
f_req_self(f_cv_Jacobian, f_cv_linear); f_req_self(f_cv_Jacobian, f_cv_linear);
f_req_children(f_cv_Jacobian, f_cvc_Jacobian); 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_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_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_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_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_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_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_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_lower_boundary);
f_req_self(f_cv_grid, f_cv_upper_boundary); f_req_self(f_cv_grid, f_cv_upper_boundary);
f_description(f_cv_lower_wall, "lower wall"); init_feature(f_cv_runave, "running average", f_type_user);
f_req_self(f_cv_lower_wall, f_cv_lower_boundary);
f_req_self(f_cv_lower_wall, f_cv_gradient);
f_description(f_cv_upper_wall, "upper wall"); init_feature(f_cv_corrfunc, "correlation function", f_type_user);
f_req_self(f_cv_upper_wall, f_cv_upper_boundary);
f_req_self(f_cv_upper_wall, f_cv_gradient);
f_description(f_cv_runave, "running average"); init_feature(f_cv_scripted, "scripted", f_type_static);
init_feature(f_cv_periodic, "periodic", f_type_static);
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");
f_req_self(f_cv_periodic, f_cv_homogeneous); f_req_self(f_cv_periodic, f_cv_homogeneous);
f_description(f_cv_scalar, "scalar"); init_feature(f_cv_scalar, "scalar", f_type_static);
f_description(f_cv_linear, "linear"); init_feature(f_cv_linear, "linear", f_type_static);
f_description(f_cv_homogeneous, "homogeneous"); init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
} }
// Initialize feature_states for each instance // Initialize feature_states for each instance
feature_states.reserve(f_cv_ntot); feature_states.reserve(f_cv_ntot);
for (i = 0; i < f_cv_ntot; i++) { 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 // Most features are available, so we set them so
// and list exceptions below // and list exceptions below
} }
// properties that may NOT be enabled as a dependency // // properties that may NOT be enabled as a dependency
int unavailable_deps[] = { // // This will be deprecated by feature types
f_cv_lower_boundary, // int unavailable_deps[] = {
f_cv_upper_boundary, // f_cv_lower_boundary,
f_cv_extended_Lagrangian, // f_cv_upper_boundary,
f_cv_Langevin, // f_cv_extended_Lagrangian,
f_cv_scripted, // f_cv_Langevin,
f_cv_periodic, // f_cv_scripted,
f_cv_scalar, // f_cv_periodic,
f_cv_linear, // f_cv_scalar,
f_cv_homogeneous // 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; // 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); 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 // The dependency below may become useful if we use dynamic atom groups
// f_req_children(f_cvc_active, f_ag_active); // 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_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_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_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 // Compute total force on first site only to avoid unwanted
// coupling to other colvars (see e.g. Ciccotti et al., 2005) // 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_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_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); 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 // 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); feature_states.reserve(f_cvc_ntot);
for (i = 0; i < colvardeps::f_cvc_ntot; i++) { 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 // Features that are implemented by all cvcs by default
// Each cvc specifies what other features are available // Each cvc specifies what other features are available
feature_states[f_cvc_active]->available = true; feature_states[f_cvc_active].available = true;
feature_states[f_cvc_gradient]->available = true; feature_states[f_cvc_gradient].available = true;
// Features that are implemented by default if their requirements are // 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 // 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_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].available = feature_states[f_cvc_scalable_com].available;
} }
@ -442,21 +460,21 @@ void colvardeps::init_ag_requires() {
features().push_back(new feature); features().push_back(new feature);
} }
f_description(f_ag_active, "active"); init_feature(f_ag_active, "active", f_type_dynamic);
f_description(f_ag_center, "translational fit"); init_feature(f_ag_center, "translational fit", f_type_static);
f_description(f_ag_rotate, "rotational fit"); init_feature(f_ag_rotate, "rotational fit", f_type_static);
f_description(f_ag_fitting_group, "reference positions group"); init_feature(f_ag_fitting_group, "reference positions group", f_type_static);
f_description(f_ag_fit_gradient_group, "fit gradient for main group"); init_feature(f_ag_fit_gradient_group, "fit gradient for main group", f_type_static);
f_description(f_ag_fit_gradient_ref, "fit gradient for reference group"); init_feature(f_ag_fit_gradient_ref, "fit gradient for reference group", f_type_static);
f_description(f_ag_atom_forces, "atomic forces"); init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
// parallel calculation implies that we have at least a scalable center of mass, // 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 // but f_ag_scalable is kept as a separate feature to deal with future dependencies
f_description(f_ag_scalable, "scalable group calculation"); init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
f_description(f_ag_scalable_com, "scalable group center of mass calculation"); 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_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_center)
// f_req_self(f_ag_min_msd_fit, f_ag_rotate) // f_req_self(f_ag_min_msd_fit, f_ag_rotate)
// f_req_exclude(f_ag_min_msd_fit, f_ag_fitting_group) // 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 // default as unavailable, not enabled
feature_states.reserve(f_ag_ntot); feature_states.reserve(f_ag_ntot);
for (i = 0; i < colvardeps::f_ag_ntot; i++) { 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 // 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 // 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) // 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; size_t i;
cvm::log("Enabled features of " + description); cvm::log("Enabled features of " + description);
for (i = 0; i < feature_states.size(); i++) { for (i = 0; i < feature_states.size(); i++) {
if (feature_states[i]->enabled) if (feature_states[i].enabled)
cvm::log("- " + features()[i]->description); cvm::log("- " + features()[i]->description);
} }
for (i=0; i<children.size(); i++) { for (i=0; i<children.size(); i++) {

86
lib/colvars/colvardeps.h
View File

@ -13,17 +13,16 @@
#include "colvarmodule.h" #include "colvarmodule.h"
#include "colvarparse.h" #include "colvarparse.h"
/// Parent class for a member object of a bias, cv or cvc etc. containing dependencies /// \brief Parent class for a member object of a bias, cv or cvc etc. containing features and
/// (features) and handling dependency resolution /// their dependencies, and handling dependency resolution
///
// Some features like colvar::f_linear have no dependencies, require() doesn't enable anything but fails if unavailable /// There are 3 kinds of features:
// Policy: those features are unavailable at all times /// 1. Dynamic features are under the control of the dependency resolution
// Other features are under user control /// system. They may be enabled or disabled depending on dependencies.
// They are unavailable unless requested by the user, then they may be enabled subject to /// 2. User features may be enabled based on user input (they may trigger a failure upon dependency resolution, though)
// satisfied dependencies /// 3. Static features are static properties of the object, determined
/// programatically at initialization time.
// It seems important to have available default to false (for safety) and enabled to false (for efficiency) ///
class colvardeps { class colvardeps {
public: public:
@ -39,12 +38,7 @@ public:
feature_state(bool a, bool e) feature_state(bool a, bool e)
: available(a), enabled(e) {} : available(a), enabled(e) {}
/// Available means: supported, subject to dependencies as listed, /// Feature may be enabled, subject to possible dependencies
/// MAY BE ENABLED AS A RESULT OF DEPENDENCY SOLVING
/// Remains false for passive flags that are set based on other object properties,
/// eg. f_cv_linear
/// Is set to true upon user request for features that are implemented by the user
/// or under his/her direct control, e.g. f_cv_scripted or f_cv_extended_Lagrangian
bool available; bool available;
/// Currently enabled - this flag is subject to change dynamically /// Currently enabled - this flag is subject to change dynamically
/// TODO consider implications for dependency solving: anyone who disables /// TODO consider implications for dependency solving: anyone who disables
@ -53,8 +47,22 @@ public:
// bool enabledOnce; // this should trigger an update when object is evaluated // bool enabledOnce; // this should trigger an update when object is evaluated
}; };
/// List of the state of all features
std::vector<feature_state *> feature_states; private:
/// List of the states of all features
std::vector<feature_state> 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 /// Describes a feature and its dependecies
/// used in a static array within each subclass /// used in a static array within each subclass
@ -83,8 +91,18 @@ public:
// features that this feature requires in children // features that this feature requires in children
std::vector<int> requires_children; std::vector<int> 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 // Accessor to array of all features with deps, static in most derived classes
// Subclasses with dynamic dependency trees may override this // Subclasses with dynamic dependency trees may override this
// with a non-static array // with a non-static array
@ -100,9 +118,8 @@ public:
/// (useful for cvcs because their children are member objects) /// (useful for cvcs because their children are member objects)
void remove_all_children(); void remove_all_children();
private: private:
// pointers to objects this object depends on // pointers to objects this object depends on
// list should be maintained by any code that modifies the object // 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 // 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 // Checks whether given feature is enabled
// Defaults to querying f_*_active // Defaults to querying f_*_active
inline bool is_enabled(int f = f_cv_active) const { 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 // Checks whether given feature is available
// Defaults to querying f_*_active // Defaults to querying f_*_active
inline bool is_available(int f = f_cv_active) const { 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 /// Parse a keyword and enable a feature accordingly
bool get_keyval_feature(colvarparse *cvp, bool get_keyval_feature(colvarparse *cvp,
@ -147,6 +176,8 @@ public:
int feature_id, bool const &def_value, int feature_id, bool const &def_value,
colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal); 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 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 // dry_run is set to true to recursively test if a feature is available, without enabling it
// int disable(int f); // int disable(int f);
@ -165,6 +196,7 @@ public:
f_cvb_get_total_force, // requires total forces f_cvb_get_total_force, // requires total forces
f_cvb_history_dependent, // depends on simulation history f_cvb_history_dependent, // depends on simulation history
f_cvb_scalar_variables, // requires scalar colvars f_cvb_scalar_variables, // requires scalar colvars
f_cvb_calc_pmf, // whether this bias will compute a PMF
f_cvb_ntot f_cvb_ntot
}; };
@ -216,12 +248,6 @@ public:
/// be used by the biases or in analysis (needs lower and upper /// be used by the biases or in analysis (needs lower and upper
/// boundary) /// boundary)
f_cv_grid, 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 /// \brief Compute running average
f_cv_runave, f_cv_runave,
/// \brief Compute time correlation function /// \brief Compute time correlation function

302
lib/colvars/colvarmodule.cpp
View File

@ -36,6 +36,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
"variable module twice.\n"); "variable module twice.\n");
return; return;
} }
depth_s = 0;
cvm::log(cvm::line_marker); cvm::log(cvm::line_marker);
cvm::log("Initializing the collective variables module, version "+ cvm::log("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n"); cvm::to_str(COLVARS_VERSION)+".\n");
@ -71,6 +74,48 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
} }
colvarmodule * colvarmodule::main()
{
return proxy->colvars;
}
std::vector<colvar *> *colvarmodule::variables()
{
return &colvars;
}
std::vector<colvar *> *colvarmodule::variables_active()
{
return &colvars_active;
}
std::vector<colvar *> *colvarmodule::variables_active_smp()
{
return &colvars_smp;
}
std::vector<int> *colvarmodule::variables_active_smp_items()
{
return &colvars_smp_items;
}
std::vector<colvarbias *> *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) int colvarmodule::read_config_file(char const *config_filename)
{ {
cvm::log(cvm::line_marker); 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) int colvarmodule::parse_config(std::string &conf)
{ {
extra_conf.clear();
// parse global options // parse global options
if (catch_input_errors(parse_global_params(conf))) { if (catch_input_errors(parse_global_params(conf))) {
return get_error(); return get_error();
@ -140,6 +206,15 @@ int colvarmodule::parse_config(std::string &conf)
return get_error(); 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(cvm::line_marker);
cvm::log("Collective variables module (re)initialized.\n"); cvm::log("Collective variables module (re)initialized.\n");
cvm::log(cvm::line_marker); 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) int colvarmodule::parse_global_params(std::string const &conf)
{ {
colvarmodule *cvm = cvm::main();
std::string index_file_name; std::string index_file_name;
if (parse->get_keyval(conf, "indexFile", 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)) { 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()) { if (colvar_conf.size()) {
cvm::log(cvm::line_marker); cvm::log(cvm::line_marker);
cvm::increase_depth(); cvm::increase_depth();
colvars.push_back(new colvar(colvar_conf)); colvars.push_back(new colvar());
if (cvm::get_error() || if (((colvars.back())->init(colvar_conf) != COLVARS_OK) ||
((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) { ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
cvm::log("Error while constructing colvar number " + cvm::log("Error while constructing colvar number " +
cvm::to_str(colvars.size()) + " : deleting."); cvm::to_str(colvars.size()) + " : deleting.");
@ -262,8 +346,7 @@ bool colvarmodule::check_new_bias(std::string &conf, char const *key)
template <class bias_type> template <class bias_type>
int colvarmodule::parse_biases_type(std::string const &conf, int colvarmodule::parse_biases_type(std::string const &conf,
char const *keyword, char const *keyword)
size_t &bias_count)
{ {
std::string bias_conf = ""; std::string bias_conf = "";
size_t conf_saved_pos = 0; size_t conf_saved_pos = 0;
@ -277,7 +360,6 @@ int colvarmodule::parse_biases_type(std::string const &conf,
return COLVARS_ERROR; return COLVARS_ERROR;
} }
cvm::decrease_depth(); cvm::decrease_depth();
bias_count++;
} else { } else {
cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n", cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n",
INPUT_ERROR); INPUT_ERROR);
@ -295,28 +377,28 @@ int colvarmodule::parse_biases(std::string const &conf)
cvm::log("Initializing the collective variables biases.\n"); cvm::log("Initializing the collective variables biases.\n");
/// initialize ABF instances /// initialize ABF instances
parse_biases_type<colvarbias_abf>(conf, "abf", n_abf_biases); parse_biases_type<colvarbias_abf>(conf, "abf");
/// initialize adaptive linear biases /// initialize adaptive linear biases
parse_biases_type<colvarbias_alb>(conf, "ALB", n_rest_biases); parse_biases_type<colvarbias_alb>(conf, "ALB");
/// initialize harmonic restraints /// initialize harmonic restraints
parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases); parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic");
/// initialize harmonic walls restraints /// initialize harmonic walls restraints
parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls", n_rest_biases); parse_biases_type<colvarbias_restraint_harmonic_walls>(conf, "harmonicWalls");
/// initialize histograms /// initialize histograms
parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases); parse_biases_type<colvarbias_histogram>(conf, "histogram");
/// initialize histogram restraints /// initialize histogram restraints
parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint", n_rest_biases); parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint");
/// initialize linear restraints /// initialize linear restraints
parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases); parse_biases_type<colvarbias_restraint_linear>(conf, "linear");
/// initialize metadynamics instances /// initialize metadynamics instances
parse_biases_type<colvarbias_meta>(conf, "metadynamics", n_meta_biases); parse_biases_type<colvarbias_meta>(conf, "metadynamics");
if (use_scripted_forces) { if (use_scripted_forces) {
cvm::log(cvm::line_marker); 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<colvarbias *>::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<colvarbias *>::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) int colvarmodule::catch_input_errors(int result)
{ {
if (result != COLVARS_OK || get_error()) { 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) { colvarbias * colvarmodule::bias_by_name(std::string const &name) {
for (std::vector<colvarbias *>::iterator bi = biases.begin(); colvarmodule *cv = cvm::main();
bi != biases.end(); for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
bi != cv->biases.end();
bi++) { bi++) {
if ((*bi)->name == name) { if ((*bi)->name == name) {
return (*bi); return (*bi);
@ -388,8 +501,9 @@ colvarbias * colvarmodule::bias_by_name(std::string const &name) {
colvar *colvarmodule::colvar_by_name(std::string const &name) { colvar *colvarmodule::colvar_by_name(std::string const &name) {
for (std::vector<colvar *>::iterator cvi = colvars.begin(); colvarmodule *cv = cvm::main();
cvi != colvars.end(); for (std::vector<colvar *>::iterator cvi = cv->colvars.begin();
cvi != cv->colvars.end();
cvi++) { cvi++) {
if ((*cvi)->name == name) { if ((*cvi)->name == name) {
return (*cvi); return (*cvi);
@ -555,9 +669,15 @@ int colvarmodule::calc_colvars()
int error_code = COLVARS_OK; int error_code = COLVARS_OK;
std::vector<colvar *>::iterator cvi; std::vector<colvar *>::iterator cvi;
// Determine which colvars are active at this time step // Determine which colvars are active at this iteration
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { variables_active()->resize(0);
(*cvi)->feature_states[colvardeps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 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 // 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 // first, calculate how much work (currently, how many active CVCs) each colvar has
colvars_smp.resize(0); variables_active_smp()->resize(0);
colvars_smp_items.resize(0); variables_active_smp_items()->resize(0);
colvars_smp.reserve(colvars.size()); variables_active_smp()->reserve(variables_active()->size());
colvars_smp_items.reserve(colvars.size()); variables_active_smp_items()->reserve(variables_active()->size());
// set up a vector containing all components // set up a vector containing all components
cvm::increase_depth(); 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(); error_code |= (*cvi)->update_cvc_flags();
size_t num_items = (*cvi)->num_active_cvcs(); size_t num_items = (*cvi)->num_active_cvcs();
colvars_smp.reserve(colvars_smp.size() + num_items); variables_active_smp()->reserve(variables_active_smp()->size() + num_items);
colvars_smp_items.reserve(colvars_smp_items.size() + num_items); variables_active_smp_items()->reserve(variables_active_smp_items()->size() + num_items);
for (size_t icvc = 0; icvc < num_items; icvc++) { for (size_t icvc = 0; icvc < num_items; icvc++) {
colvars_smp.push_back(*cvi); variables_active_smp()->push_back(*cvi);
colvars_smp_items.push_back(icvc); variables_active_smp_items()->push_back(icvc);
} }
} }
cvm::decrease_depth(); cvm::decrease_depth();
@ -592,8 +711,7 @@ int colvarmodule::calc_colvars()
error_code |= proxy->smp_colvars_loop(); error_code |= proxy->smp_colvars_loop();
cvm::increase_depth(); 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)->collect_cvc_data(); error_code |= (*cvi)->collect_cvc_data();
} }
cvm::decrease_depth(); cvm::decrease_depth();
@ -602,8 +720,7 @@ int colvarmodule::calc_colvars()
// calculate colvars one at a time // calculate colvars one at a time
cvm::increase_depth(); 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)->calc(); error_code |= (*cvi)->calc();
if (cvm::get_error()) { if (cvm::get_error()) {
return COLVARS_ERROR; return COLVARS_ERROR;
@ -627,6 +744,18 @@ int colvarmodule::calc_biases()
std::vector<colvarbias *>::iterator bi; std::vector<colvarbias *>::iterator bi;
int error_code = COLVARS_OK; 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 SMP support is available, split up the work
if (proxy->smp_enabled() == COLVARS_OK) { if (proxy->smp_enabled() == COLVARS_OK) {
@ -645,7 +774,7 @@ int colvarmodule::calc_biases()
} }
cvm::increase_depth(); 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(); error_code |= (*bi)->update();
if (cvm::get_error()) { if (cvm::get_error()) {
return error_code; return error_code;
@ -654,12 +783,10 @@ int colvarmodule::calc_biases()
cvm::decrease_depth(); cvm::decrease_depth();
} }
cvm::real total_bias_energy = 0.0; for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
for (bi = biases.begin(); bi != biases.end(); bi++) {
total_bias_energy += (*bi)->get_energy(); total_bias_energy += (*bi)->get_energy();
} }
proxy->add_energy(total_bias_energy);
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
} }
@ -675,7 +802,7 @@ int colvarmodule::update_colvar_forces()
if (cvm::debug() && biases.size()) if (cvm::debug() && biases.size())
cvm::log("Collecting forces from all biases.\n"); cvm::log("Collecting forces from all biases.\n");
cvm::increase_depth(); cvm::increase_depth();
for (bi = biases.begin(); bi != biases.end(); bi++) { for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
(*bi)->communicate_forces(); (*bi)->communicate_forces();
if (cvm::get_error()) { if (cvm::get_error()) {
return COLVARS_ERROR; return COLVARS_ERROR;
@ -687,6 +814,11 @@ int colvarmodule::update_colvar_forces()
error_code |= calc_scripted_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; cvm::real total_colvar_energy = 0.0;
// sum up the forces for each colvar, including wall forces // sum up the forces for each colvar, including wall forces
// and integrate any internal // and integrate any internal
@ -695,7 +827,7 @@ int colvarmodule::update_colvar_forces()
cvm::log("Updating the internal degrees of freedom " cvm::log("Updating the internal degrees of freedom "
"of colvars (if they have any).\n"); "of colvars (if they have any).\n");
cvm::increase_depth(); 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 // Here we call even inactive colvars, so they accumulate biasing forces
// as well as update their extended-system dynamics // as well as update their extended-system dynamics
total_colvar_energy += (*cvi)->update_forces_energy(); total_colvar_energy += (*cvi)->update_forces_energy();
@ -704,6 +836,8 @@ int colvarmodule::update_colvar_forces()
} }
} }
cvm::decrease_depth(); 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); proxy->add_energy(total_colvar_energy);
// make collective variables communicate their forces to their // make collective variables communicate their forces to their
@ -711,9 +845,8 @@ int colvarmodule::update_colvar_forces()
if (cvm::debug()) if (cvm::debug())
cvm::log("Communicating forces from the colvars to the atoms.\n"); cvm::log("Communicating forces from the colvars to the atoms.\n");
cvm::increase_depth(); 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(colvardeps::f_cv_gradient)) {
if (!(*cvi)->is_enabled()) continue;
(*cvi)->communicate_forces(); (*cvi)->communicate_forces();
if (cvm::get_error()) { if (cvm::get_error()) {
return COLVARS_ERROR; return COLVARS_ERROR;
@ -800,8 +933,8 @@ int colvarmodule::analyze()
cvm::log("Performing analysis.\n"); cvm::log("Performing analysis.\n");
// perform colvar-specific analysis // perform colvar-specific analysis
for (std::vector<colvar *>::iterator cvi = colvars.begin(); for (std::vector<colvar *>::iterator cvi = variables_active()->begin();
cvi != colvars.end(); cvi != variables_active()->end();
cvi++) { cvi++) {
cvm::increase_depth(); cvm::increase_depth();
(*cvi)->analyze(); (*cvi)->analyze();
@ -825,8 +958,8 @@ int colvarmodule::setup()
{ {
if (this->size() == 0) return cvm::get_error(); if (this->size() == 0) return cvm::get_error();
// loop over all components of all colvars to reset masses of all groups // loop over all components of all colvars to reset masses of all groups
for (std::vector<colvar *>::iterator cvi = colvars.begin(); for (std::vector<colvar *>::iterator cvi = variables()->begin();
cvi != colvars.end(); cvi++) { cvi != variables()->end(); cvi++) {
(*cvi)->setup(); (*cvi)->setup();
} }
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); 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"); cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
} }
output_prefix = proxy->output_prefix(); output_prefix() = proxy->output_prefix();
if (output_prefix.size()) { if (output_prefix().size()) {
cvm::log("The final output state file will be \""+ cvm::log("The final output state file will be \""+
(output_prefix.size() ? (output_prefix().size() ?
std::string(output_prefix+".colvars.state") : std::string(output_prefix()+".colvars.state") :
std::string("colvars.state"))+"\".\n"); std::string("colvars.state"))+"\".\n");
// cvm::log (cvm::line_marker); // cvm::log (cvm::line_marker);
} }
cv_traj_name = cv_traj_name =
(output_prefix.size() ? (output_prefix().size() ?
std::string(output_prefix+".colvars.traj") : std::string(output_prefix()+".colvars.traj") :
std::string("")); std::string(""));
if (cv_traj_freq && cv_traj_name.size()) { if (cv_traj_freq && cv_traj_name.size()) {
@ -1077,14 +1210,14 @@ continue the previous simulation.\n\n");
cvm::log(cvm::line_marker); cvm::log(cvm::line_marker);
// update this ahead of time in this special case // update this ahead of time in this special case
output_prefix = proxy->input_prefix(); output_prefix() = proxy->input_prefix();
cvm::log("All output files will now be saved with the prefix \""+output_prefix+".tmp.*\".\n"); cvm::log("All output files will now be saved with the prefix \""+output_prefix()+".tmp.*\".\n");
cvm::log(cvm::line_marker); cvm::log(cvm::line_marker);
cvm::log("Please review the important warning above. After that, you may rename:\n\ 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\ to:\n\
\""+ proxy->input_prefix()+".colvars.state\"\n"); \""+ proxy->input_prefix()+".colvars.state\"\n");
output_prefix = output_prefix+".tmp"; output_prefix() = output_prefix()+".tmp";
write_output_files(); write_output_files();
cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR); 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 // if this is a simulation run (i.e. not a postprocessing), output data
// must be written to be able to restart the simulation // must be written to be able to restart the simulation
std::string const out_name = std::string const out_name =
(output_prefix.size() ? (output_prefix().size() ?
std::string(output_prefix+".colvars.state") : std::string(output_prefix()+".colvars.state") :
std::string("colvars.state")); std::string("colvars.state"));
cvm::log("Saving collective variables state to \""+out_name+"\".\n"); 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) 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) if (d > 0)
proxy->log((std::string(2*d, ' '))+message); proxy->log((std::string(2*d, ' '))+message);
else else
@ -1365,19 +1499,20 @@ void cvm::decrease_depth()
size_t & cvm::depth() size_t & cvm::depth()
{ {
// NOTE: do not call log() or error() here, to avoid recursion // 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 (proxy->smp_enabled() == COLVARS_OK) {
if (depth_v.size() != nt) { int const nt = proxy->smp_num_threads();
// update array of depths if (int(cv->depth_v.size()) != nt) {
proxy->smp_lock(); proxy->smp_lock();
if (depth_v.size() > 0) { depth_s = depth_v[0]; } // update array of depths
depth_v.clear(); if (cv->depth_v.size() > 0) { cv->depth_s = cv->depth_v[0]; }
depth_v.assign(nt, depth_s); cv->depth_v.clear();
cv->depth_v.assign(nt, cv->depth_s);
proxy->smp_unlock(); 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); FILE_ERROR);
} }
} }
cvm::index_group_names.push_back(group_name); index_group_names.push_back(group_name);
cvm::index_groups.push_back(std::vector<int> ()); index_groups.push_back(std::vector<int>());
} else { } else {
cvm::error("Error: in parsing index file \""+ cvm::error("Error: in parsing index file \""+
std::string(filename)+"\".\n", std::string(filename)+"\".\n",
@ -1462,7 +1597,7 @@ int cvm::read_index_file(char const *filename)
int atom_number = 1; int atom_number = 1;
size_t pos = is.tellg(); size_t pos = is.tellg();
while ( (is >> atom_number) && (atom_number > 0) ) { 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(); pos = is.tellg();
} }
is.clear(); is.clear();
@ -1532,8 +1667,8 @@ int cvm::load_coords_xyz(char const *filename,
+ std::string(filename) + ".\n", INPUT_ERROR); + std::string(filename) + ".\n", INPUT_ERROR);
} }
// skip comment line // skip comment line
std::getline(xyz_is, line); cvm::getline(xyz_is, line);
std::getline(xyz_is, line); cvm::getline(xyz_is, line);
xyz_is.width(255); xyz_is.width(255);
std::vector<atom_pos>::iterator pos_i = pos.begin(); std::vector<atom_pos>::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++) { for ( ; pos_i != pos.end() ; pos_i++, index++) {
while ( next < *index ) { while ( next < *index ) {
std::getline(xyz_is, line); cvm::getline(xyz_is, line);
next++; next++;
} }
xyz_is >> symbol; xyz_is >> symbol;
@ -1558,15 +1693,9 @@ int cvm::load_coords_xyz(char const *filename,
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
} }
// static pointers
std::vector<colvar *> colvarmodule::colvars;
std::vector<colvarbias *> 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 // static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07; cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
@ -1575,14 +1704,9 @@ long colvarmodule::it = 0;
long colvarmodule::it_restart = 0; long colvarmodule::it_restart = 0;
size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::restart_out_freq = 0;
size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::cv_traj_freq = 0;
size_t colvarmodule::depth_s = 0;
std::vector<size_t> colvarmodule::depth_v(0);
bool colvarmodule::b_analysis = false; bool colvarmodule::b_analysis = false;
std::list<std::string> colvarmodule::index_group_names; bool colvarmodule::use_scripted_forces = false;
std::list<std::vector<int> > colvarmodule::index_groups; bool colvarmodule::scripting_after_biases = true;
bool colvarmodule::use_scripted_forces = false;
bool colvarmodule::scripting_after_biases = true;
std::string colvarmodule::output_prefix = "";
// i/o constants // i/o constants
size_t const colvarmodule::it_width = 12; size_t const colvarmodule::it_width = 12;

115
lib/colvars/colvarmodule.h
View File

@ -11,7 +11,7 @@
#define COLVARMODULE_H #define COLVARMODULE_H
#ifndef COLVARS_VERSION #ifndef COLVARS_VERSION
#define COLVARS_VERSION "2016-12-27" #define COLVARS_VERSION "2017-03-09"
#endif #endif
#ifndef COLVARS_DEBUG #ifndef COLVARS_DEBUG
@ -161,12 +161,37 @@ public:
/// dt) /// dt)
static real debug_gradients_step_size; static real debug_gradients_step_size;
/// Prefix for all output files for this run private:
static std::string output_prefix;
/// 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 /// Array of collective variables
static std::vector<colvar *> colvars; std::vector<colvar *> colvars;
/// Array of collective variables
std::vector<colvar *> colvars_active;
/// Collective variables to be calculated on different threads;
/// colvars with multple items (e.g. multiple active CVCs) are duplicated
std::vector<colvar *> colvars_smp;
/// Indexes of the items to calculate for each colvar
std::vector<int> colvars_smp_items;
public:
/// Array of collective variables
std::vector<colvar *> *variables();
/* TODO: implement named CVCs /* TODO: implement named CVCs
/// Array of named (reusable) collective variable components /// Array of named (reusable) collective variable components
@ -177,26 +202,31 @@ public:
} }
*/ */
/// Collective variables with the active flag on
std::vector<colvar *> *variables_active();
/// Collective variables to be calculated on different threads; /// Collective variables to be calculated on different threads;
/// colvars with multple items (e.g. multiple active CVCs) are duplicated /// colvars with multple items (e.g. multiple active CVCs) are duplicated
std::vector<colvar *> colvars_smp; std::vector<colvar *> *variables_active_smp();
/// Indexes of the items to calculate for each colvar /// Indexes of the items to calculate for each colvar
std::vector<int> colvars_smp_items; std::vector<int> *variables_active_smp_items();
/// Array of collective variable biases /// Array of collective variable biases
static std::vector<colvarbias *> biases; std::vector<colvarbias *> biases;
/// \brief Number of ABF biases initialized (in normal conditions
/// should be 1) /// Energy of built-in and scripted biases, summed per time-step
static size_t n_abf_biases; real total_bias_energy;
/// \brief Number of metadynamics biases initialized (in normal
/// conditions should be 1) private:
static size_t n_meta_biases;
/// \brief Number of restraint biases initialized (no limit on the /// Array of active collective variable biases
/// number) std::vector<colvarbias *> biases_active_;
static size_t n_rest_biases;
/// \brief Number of histograms initialized (no limit on the public:
/// number)
static size_t n_histo_biases; /// Array of active collective variable biases
std::vector<colvarbias *> *biases_active();
/// \brief Whether debug output should be enabled (compile-time option) /// \brief Whether debug output should be enabled (compile-time option)
static inline bool debug() static inline bool debug()
@ -205,10 +235,7 @@ public:
} }
/// \brief How many objects are configured yet? /// \brief How many objects are configured yet?
inline size_t size() const size_t size() const;
{
return colvars.size() + biases.size();
}
/// \brief Constructor \param config_name Configuration file name /// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name /// \param restart_name (optional) Restart file name
@ -230,9 +257,12 @@ public:
/// \brief Parse a "clean" config string (no comments) /// \brief Parse a "clean" config string (no comments)
int parse_config(std::string &conf); int parse_config(std::string &conf);
// Parse functions (setup internal data based on a string) // 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 /// Parse the few module's global parameters
int parse_global_params(std::string const &conf); int parse_global_params(std::string const &conf);
@ -242,14 +272,33 @@ public:
/// Parse and initialize collective variable biases /// Parse and initialize collective variable biases
int parse_biases(std::string const &conf); 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 /// Parse and initialize collective variable biases of a specific type
template <class bias_type> template <class bias_type>
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 /// Test error condition and keyword parsing
/// on error, delete new bias /// on error, delete new bias
bool check_new_bias(std::string &conf, char const *key); 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: private:
/// Useful wrapper to interrupt parsing if any error occurs /// Useful wrapper to interrupt parsing if any error occurs
int catch_input_errors(int result); 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 /// \brief Names of groups from a Gromacs .ndx file to be read at startup
static std::list<std::string> index_group_names; std::list<std::string> index_group_names;
/// \brief Groups from a Gromacs .ndx file read at startup /// \brief Groups from a Gromacs .ndx file read at startup
static std::list<std::vector<int> > index_groups; std::list<std::vector<int> > index_groups;
/// \brief Read a Gromacs .ndx file /// \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 /// \brief Create atoms from a file \param filename name of the file
@ -515,13 +564,13 @@ protected:
/// Output restart file /// Output restart file
colvarmodule::ofstream restart_out_os; colvarmodule::ofstream restart_out_os;
protected: private:
/// Counter for the current depth in the object hierarchy (useg e.g. in output) /// 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 /// Thread-specific depth
static std::vector<size_t> depth_v; std::vector<size_t> depth_v;
public: public:
@ -552,6 +601,10 @@ public:
/// from the hosting program; it is static in order to be accessible /// from the hosting program; it is static in order to be accessible
/// from static functions in the colvarmodule class /// from static functions in the colvarmodule class
static colvarproxy *proxy; static colvarproxy *proxy;
/// \brief Accessor for the above
static colvarmodule *main();
}; };

7
lib/colvars/colvarparse.cpp
View File

@ -503,7 +503,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
std::string line; std::string line;
std::istringstream is(conf); std::istringstream is(conf);
while (std::getline(is, line)) { while (cvm::getline(is, line)) {
if (line.size() == 0) if (line.size() == 0)
continue; continue;
if (line.find_first_not_of(white_space) == 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::istream & colvarparse::getline_nocomments(std::istream &is,
std::string &line, std::string &line)
char const delim)
{ {
std::getline(is, line, delim); cvm::getline(is, line);
size_t const comment = line.find('#'); size_t const comment = line.find('#');
if (comment != std::string::npos) { if (comment != std::string::npos) {
line.erase(comment); line.erase(comment);

3
lib/colvars/colvarparse.h
View File

@ -304,8 +304,7 @@ public:
/// \brief Works as std::getline() but also removes everything /// \brief Works as std::getline() but also removes everything
/// between a comment character and the following newline /// between a comment character and the following newline
static std::istream & getline_nocomments(std::istream &is, static std::istream & getline_nocomments(std::istream &is,
std::string &s, std::string &s);
char const delim = '\n');
/// Check if the content of the file has matching braces /// Check if the content of the file has matching braces
bool brace_check(std::string const &conf, bool brace_check(std::string const &conf,

2
lib/colvars/colvarproxy.h
View File

@ -336,8 +336,6 @@ public:
/// Are total forces being used? /// Are total forces being used?
virtual bool total_forces_enabled() const virtual bool total_forces_enabled() const
{ {
cvm::error("Error: total forces are currently not implemented.\n",
COLVARS_NOT_IMPLEMENTED);
return false; return false;
} }

117
lib/colvars/colvarscript.cpp
View File

@ -12,6 +12,7 @@
#include <string.h> #include <string.h>
#include "colvarscript.h" #include "colvarscript.h"
#include "colvardeps.h"
colvarscript::colvarscript(colvarproxy *p) 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(); result = "Syntax error\n" + help_string();
return COLVARSCRIPT_ERROR; return COLVARSCRIPT_ERROR;
} }
@ -263,10 +274,11 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
} }
if (subcmd == "delete") { if (subcmd == "delete") {
if (cv->biases.size() > 0) { size_t i;
result = "Cannot delete a colvar currently used by biases, delete those biases first"; for (i = 0; i < cv->biases.size(); i++) {
return COLVARSCRIPT_ERROR; delete cv->biases[i];
} }
cv->biases.resize(0);
// colvar destructor is tasked with the cleanup // colvar destructor is tasked with the cleanup
delete cv; delete cv;
// TODO this could be done by the destructors // TODO this could be done by the destructors
@ -338,9 +350,8 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
return COLVARS_OK; return COLVARS_OK;
} }
if (subcmd == "state") { if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
cv->print_state(); return proc_features(cv, argc, argv);
return COLVARS_OK;
} }
result = "Syntax error\n" + help_string(); result = "Syntax error\n" + help_string();
@ -379,11 +390,6 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
return COLVARS_OK; return COLVARS_OK;
} }
if (subcmd == "state") {
b->print_state();
return COLVARS_OK;
}
// Subcommands for MW ABF // Subcommands for MW ABF
if (subcmd == "bin") { if (subcmd == "bin") {
int r = b->current_bin(); int r = b->current_bin();
@ -420,6 +426,10 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
return COLVARS_OK; return COLVARS_OK;
} }
if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
return proc_features(b, argc, argv);
}
if (argc >= 4) { if (argc >= 4) {
std::string param = argv[3]; std::string param = argv[3];
if (subcmd == "count") { 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<colvardeps::feature *> &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 colvarscript::help_string()
{ {
std::string buf; std::string buf;
@ -459,6 +546,7 @@ Input and output:\n\
load <file name> -- load a state file (requires configuration)\n\ load <file name> -- load a state file (requires configuration)\n\
save <file name> -- save a state file (requires configuration)\n\ save <file name> -- save a state file (requires configuration)\n\
update -- recalculate colvars and biases\n\ update -- recalculate colvars and biases\n\
addenergy <E> -- add <E> to the total bias energy\n\
printframe -- return a summary of the current frame\n\ printframe -- return a summary of the current frame\n\
printframelabels -- return labels to annotate printframe's output\n"; printframelabels -- return labels to annotate printframe's output\n";
@ -478,12 +566,17 @@ Accessing collective variables:\n\
colvar <name> addforce <F> -- apply given force on colvar <name>\n\ colvar <name> addforce <F> -- apply given force on colvar <name>\n\
colvar <name> getconfig -- return config string of colvar <name>\n\ colvar <name> getconfig -- return config string of colvar <name>\n\
colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\ colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\
colvar <name> get <f> -- get the value of the colvar feature <f>\n\
colvar <name> set <f> <val> -- set the value of the colvar feature <f>\n\
\n\ \n\
Accessing biases:\n\ Accessing biases:\n\
bias <name> energy -- return the current energy of bias <name>\n\ bias <name> energy -- return the current energy of bias <name>\n\
bias <name> update -- recalculate bias <name>\n\ bias <name> update -- recalculate bias <name>\n\
bias <name> delete -- delete bias <name>\n\ bias <name> delete -- delete bias <name>\n\
bias <name> getconfig -- return config string of bias <name>\n"; bias <name> getconfig -- return config string of bias <name>\n\
bias <name> get <f> -- get the value of the bias feature <f>\n\
bias <name> set <f> <val> -- set the value of the bias feature <f>\n\
";
return buf; return buf;
} }

4
lib/colvars/colvarscript.h
View File

@ -51,6 +51,10 @@ private:
/// Run subcommands on bias /// Run subcommands on bias
int proc_bias(int argc, char const *argv[]); 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 /// Builds and return a short help
std::string help_string(void); std::string help_string(void);
}; };

3
src/KOKKOS/comm_kokkos.cpp
View File

@ -93,6 +93,9 @@ CommKokkos::~CommKokkos()
void CommKokkos::init() void CommKokkos::init()
{ {
maxsend = BUFMIN;
maxrecv = BUFMIN;
grow_send_kokkos(maxsend+bufextra,0,Host); grow_send_kokkos(maxsend+bufextra,0,Host);
grow_recv_kokkos(maxrecv,Host); grow_recv_kokkos(maxrecv,Host);

19
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -352,12 +352,17 @@ int colvarproxy_lammps::smp_enabled()
int colvarproxy_lammps::smp_colvars_loop() int colvarproxy_lammps::smp_colvars_loop()
{ {
colvarmodule *cv = this->colvars; colvarmodule *cv = this->colvars;
colvarproxy_lammps *proxy = (colvarproxy_lammps *) cv->proxy;
#pragma omp parallel for #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()) { 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(); return cvm::get_error();
} }
@ -367,11 +372,13 @@ int colvarproxy_lammps::smp_biases_loop()
{ {
colvarmodule *cv = this->colvars; colvarmodule *cv = this->colvars;
#pragma omp parallel for #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()) { 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(); return cvm::get_error();
} }

2
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -29,7 +29,7 @@
#endif #endif
#ifndef COLVARPROXY_VERSION #ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2016-12-27" #define COLVARPROXY_VERSION "2017-01-09"
#endif #endif
/* struct for packed data communication of coordinates and forces. */ /* struct for packed data communication of coordinates and forces. */

21
src/balance.cpp
View File

@ -456,6 +456,7 @@ void Balance::options(int iarg, int narg, char **arg)
wtflag = 0; wtflag = 0;
varflag = 0; varflag = 0;
oldrcb = 0;
outflag = 0; outflag = 0;
int outarg = 0; int outarg = 0;
fp = NULL; fp = NULL;
@ -491,6 +492,9 @@ void Balance::options(int iarg, int narg, char **arg)
} }
iarg += 2+nopt; iarg += 2+nopt;
} else if (strcmp(arg[iarg],"old") == 0) {
oldrcb = 1;
iarg++;
} else if (strcmp(arg[iarg],"out") == 0) { } else if (strcmp(arg[iarg],"out") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal (fix) balance command"); if (iarg+2 > narg) error->all(FLERR,"Illegal (fix) balance command");
outflag = 1; outflag = 1;
@ -641,11 +645,20 @@ int *Balance::bisection(int sortflag)
// invoke RCB // invoke RCB
// then invert() to create list of proc assignments for my atoms // 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) { if (oldrcb) {
weight = fixstore->vstore; if (wtflag) {
rcb->compute(dim,atom->nlocal,atom->x,weight,shrinklo,shrinkhi); weight = fixstore->vstore;
} else rcb->compute(dim,atom->nlocal,atom->x,NULL,shrinklo,shrinkhi); 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); rcb->invert(sortflag);

1
src/balance.h
View File

@ -53,6 +53,7 @@ class Balance : protected Pointers {
int style; // style of LB int style; // style of LB
int xflag,yflag,zflag; // xyz LB flags int xflag,yflag,zflag; // xyz LB flags
double *user_xsplit,*user_ysplit,*user_zsplit; // params for xyz LB double *user_xsplit,*user_ysplit,*user_zsplit; // params for xyz LB
int oldrcb; // use old-style RCB compute
int nitermax; // params for shift LB int nitermax; // params for shift LB
double stopthresh; double stopthresh;

57
src/fix_halt.cpp
View File

@ -30,9 +30,10 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{BONDMAX,VARIABLE}; enum{BONDMAX,TLIMIT,VARIABLE};
enum{LT,LE,GT,GE,EQ,NEQ,XOR}; enum{LT,LE,GT,GE,EQ,NEQ,XOR};
enum{HARD,SOFT,CONTINUE}; enum{HARD,SOFT,CONTINUE};
enum{NOMSG,YESMSG};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -47,7 +48,8 @@ FixHalt::FixHalt(LAMMPS *lmp, int narg, char **arg) :
idvar = NULL; 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) { else if (strncmp(arg[4],"v_",2) == 0) {
attribute = VARIABLE; attribute = VARIABLE;
int n = strlen(arg[4]); int n = strlen(arg[4]);
@ -73,6 +75,7 @@ FixHalt::FixHalt(LAMMPS *lmp, int narg, char **arg) :
// parse optional args // parse optional args
eflag = SOFT; eflag = SOFT;
msgflag = YESMSG;
int iarg = 7; int iarg = 7;
while (iarg < narg) { 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 if (strcmp(arg[iarg+1],"continue") == 0) eflag = CONTINUE;
else error->all(FLERR,"Illegal fix halt command"); else error->all(FLERR,"Illegal fix halt command");
iarg += 2; 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"); } else error->all(FLERR,"Illegal fix halt command");
} }
@ -125,6 +134,11 @@ void FixHalt::init()
if (input->variable->equalstyle(ivar) == 0) if (input->variable->equalstyle(ivar) == 0)
error->all(FLERR,"Fix halt variable is not equal-style variable"); 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; double attvalue;
if (attribute == BONDMAX) attvalue = bondmax(); if (attribute == TLIMIT) {
else { if (update->ntimestep != nextstep) return;
attvalue = tlimit();
} else if (attribute == BONDMAX) {
attvalue = bondmax();
} else {
modify->clearstep_compute(); modify->clearstep_compute();
attvalue = input->variable->compute_equal(ivar); attvalue = input->variable->compute_equal(ivar);
modify->addstep_compute(update->ntimestep + nevery); 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", sprintf(str,"Fix halt %s condition met on step %ld with value %g",
id,update->ntimestep,attvalue); id,update->ntimestep,attvalue);
if (eflag == HARD) error->all(FLERR,str); if (eflag == HARD) {
else if (eflag == SOFT || eflag == CONTINUE) { error->all(FLERR,str);
if (comm->me == 0) error->message(FLERR,str); } else if (eflag == SOFT || eflag == CONTINUE) {
if (comm->me == 0 && msgflag == YESMSG) error->message(FLERR,str);
timer->force_timeout(); timer->force_timeout();
} }
} }
@ -218,3 +237,27 @@ double FixHalt::bondmax()
return sqrt(maxall); 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<bigint> (tratio*value/cpu * elapsed);
nextstep = (final/nevery)*nevery + nevery;
tratio = 1.0;
}
//printf("EVAL %ld %g %d\n",update->ntimestep,cpu,nevery);
return cpu;
}

6
src/fix_halt.h
View File

@ -35,11 +35,13 @@ class FixHalt : public Fix {
void post_run(); void post_run();
private: private:
int attribute,operation,eflag,ivar; int attribute,operation,eflag,msgflag,ivar;
double value; bigint nextstep;
double value,tratio;
char *idvar; char *idvar;
double bondmax(); double bondmax();
double tlimit();
}; };
} }

7
src/neighbor.cpp
View File

@ -1129,6 +1129,11 @@ void Neighbor::morph_halffull()
if (!irq->half) continue; 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 // these lists are created other ways, no need for halffull
// do want to process skip lists // do want to process skip lists
@ -1154,8 +1159,6 @@ void Neighbor::morph_halffull()
// this includes custom cutoff set by requestor // this includes custom cutoff set by requestor
// no need to check respaouter b/c it stores same pairs // no need to check respaouter b/c it stores same pairs
// no need to check dnum b/c only set for history // 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->ghost != jrq->ghost) continue;
if (irq->size != jrq->size) continue; if (irq->size != jrq->size) continue;

523
src/rcb.cpp
View File

@ -42,7 +42,7 @@ RCB::RCB(LAMMPS *lmp) : Pointers(lmp)
dots = NULL; dots = NULL;
nlist = maxlist = 0; nlist = maxlist = 0;
dotlist = dotmark = NULL; dotlist = dotmark = dotmark_select = NULL;
maxbuf = 0; maxbuf = 0;
buf = NULL; buf = NULL;
@ -73,6 +73,7 @@ RCB::~RCB()
memory->sfree(dots); memory->sfree(dots);
memory->destroy(dotlist); memory->destroy(dotlist);
memory->destroy(dotmark); memory->destroy(dotmark);
memory->destroy(dotmark_select);
memory->sfree(buf); memory->sfree(buf);
memory->destroy(recvproc); memory->destroy(recvproc);
@ -91,6 +92,9 @@ RCB::~RCB()
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
perform RCB balancing of N particles at coords X in bounding box LO/HI 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 = NULL, ignore per-particle weights
if wt defined, per-particle weights > 0.0 if wt defined, per-particle weights > 0.0
dimension = 2 or 3 dimension = 2 or 3
@ -103,6 +107,523 @@ RCB::~RCB()
void RCB::compute(int dimension, int n, double **x, double *wt, void RCB::compute(int dimension, int n, double **x, double *wt,
double *bboxlo, double *bboxhi) 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 i,j,k;
int keep,outgoing,incoming,incoming2; int keep,outgoing,incoming,incoming2;

2
src/rcb.h
View File

@ -42,6 +42,7 @@ class RCB : protected Pointers {
RCB(class LAMMPS *); RCB(class LAMMPS *);
~RCB(); ~RCB();
void compute(int, int, double **, double *, double *, double *); void compute(int, int, double **, double *, double *, double *);
void compute_old(int, int, double **, double *, double *, double *);
void invert(int sortflag = 0); void invert(int sortflag = 0);
bigint memory_usage(); bigint memory_usage();
@ -99,6 +100,7 @@ class RCB : protected Pointers {
int maxlist; int maxlist;
int *dotlist; int *dotlist;
int *dotmark; int *dotmark;
int *dotmark_select;
int maxbuf; int maxbuf;
Dot *buf; Dot *buf;

8
src/thermo.cpp
View File

@ -55,7 +55,8 @@ using namespace MathConst;
// customize a new keyword by adding to this list: // 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 // atoms, temp, press, pe, ke, etotal, enthalpy
// evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail // evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail
// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz, // 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); 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;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "7 Mar 2017" #define LAMMPS_VERSION "10 Mar 2017"