Merge pull request #408 from giacomofiorin/colvars-update-2017-03-10 (470353e3) · Commits · 郑智淋 / lammps

doc/src/PDF/colvars-refman-lammps.pdf

+677 B (553 KiB)

File changed.

No diff preview for this file type.

lib/colvars/colvar.cpp

+191 −171

Original line number	Original line	Diff line number	Diff line
	// -- 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;

			colvarmodule *cv = cvm::main();

	get_keyval(conf, "name", this->name,		get_keyval(conf, "name", this->name,
	(std::string("colvar")+cvm::to_str(cvm::colvars.size()+1)));		(std::string("colvar")+cvm::to_str(cv->variables()->size()+1)));

	if (cvm::colvar_by_name(this->name) != NULL) {		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_kwidthwidth)+"\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_kwidthwidth)+"\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,16 +356,21 @@ 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;
			}

			return COLVARS_OK;
	}		}

	get_keyval(conf, "timeStepFactor", time_step_factor, 1);

			int colvar::init_extended_Lagrangian(std::string const &conf)
	{		{
	bool b_extended_Lagrangian;		bool b_extended_Lagrangian;
	get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);		get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
	@@ -309,13 +381,7 @@ colvar::colvar(std::string const &conf)
	cvm::log("Enabling the extended Lagrangian term for colvar \""+		cvm::log("Enabling the extended Lagrangian term for colvar \""+
	this->name+"\".\n");		this->name+"\".\n");

	// Make feature available only on user request
	provide(f_cv_extended_Lagrangian);
	enable(f_cv_extended_Lagrangian);		enable(f_cv_extended_Lagrangian);
	provide(f_cv_Langevin);

	// The extended mass will apply forces
	enable(f_cv_gradient);

	xr.type(value());		xr.type(value());
	vr.type(value());		vr.type(value());
	@@ -329,11 +395,13 @@ colvar::colvar(std::string const &conf)
	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);		get_keyval(conf, "extendedFluctuation", tolerance);
	if (tolerance <= 0.0) {		if (tolerance <= 0.0) {
	cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);		cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR);
			return INPUT_ERROR;
	}		}
	ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);		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");		cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
	@@ -357,6 +425,7 @@ colvar::colvar(std::string const &conf)
	get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);		get_keyval(conf, "extendedLangevinDamping", ext_gamma, 1.0);
	if (ext_gamma < 0.0) {		if (ext_gamma < 0.0) {
	cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);		cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR);
			return INPUT_ERROR;
	}		}
	if (ext_gamma != 0.0) {		if (ext_gamma != 0.0) {
	enable(f_cv_Langevin);		enable(f_cv_Langevin);
	@@ -364,8 +433,13 @@ colvar::colvar(std::string const &conf)
	ext_sigma = std::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());		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;

lib/colvars/colvar.h

+13 −1

Original line number	Original line	Diff line number	Diff line
	@@ -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,

lib/colvars/colvaratoms.cpp

+8 −6

Original line number	Original line	Diff line number	Diff line
	@@ -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);

lib/colvars/colvarbias.cpp

+26 −25

Original line number	Original line	Diff line number	Diff line
	@@ -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]);
	}		}
	}		}

Admin message