colvar.C
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Nov 18, 06:56

colvar.C
View Options

	// -- c++ --

	#include "colvarmodule.h"
	#include "colvarvalue.h"
	#include "colvarparse.h"
	#include "colvar.h"
	#include "colvarcomp.h"
	#include <algorithm>


	// XX TODO make the acf code handle forces as well as values and velocities


	colvar::colvar (std::string const &conf)
	{
	cvm::log ("Initializing a new collective variable.\n");

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

	for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
	cvi < cvm::colvars.end();
	cvi++) {
	if ((*cvi)->name == this->name)
	cvm::fatal_error ("Error: this colvar has the same name, \""+this->name+
	"\", as another colvar.\n");
	}

	// all tasks disabled by default
	for (size_t i = 0; i < task_ntot; i++) {
	tasks[i] = false;
	}

	kinetic_energy = 0.0;
	potential_energy = 0.0;

	// read the configuration and set up corresponding instances, for
	// each type of component implemented
	#define initialize_components(def_desc,def_config_key,def_class_name) \
	{ \
	size_t def_count = 0; \
	std::string def_conf = ""; \
	size_t pos = 0; \
	while ( this->key_lookup (conf, \
	def_config_key, \
	def_conf, \
	pos) ) { \
	if (!def_conf.size()) continue; \
	cvm::log ("Initializing " \
	"a new \""+std::string (def_config_key)+"\" component"+ \
	(cvm::debug() ? ", with configuration:\n"+def_conf \
	: ".\n")); \
	cvm::increase_depth(); \
	cvc *cvcp = new colvar::def_class_name (def_conf); \
	if (cvcp != NULL) { \
	cvcs.push_back (cvcp); \
	cvcp->check_keywords (def_conf, def_config_key); \
	cvm::decrease_depth(); \
	} else { \
	cvm::fatal_error ("Error: in allocating component \"" \
	def_config_key"\".\n"); \
	} \
	if ( (cvcp->period != 0.0) \|\| (cvcp->wrap_center != 0.0) ) { \
	if ( (cvcp->function_type != std::string ("distance_z")) && \
	(cvcp->function_type != std::string ("dihedral")) && \
	(cvcp->function_type != std::string ("spin_angle")) ) { \
	cvm::fatal_error ("Error: invalid use of period and/or " \
	"wrapAround in a \""+ \
	std::string (def_config_key)+ \
	"\" component.\n"+ \
	"Period: "+cvm::to_str(cvcp->period) + \
	" wrapAround: "+cvm::to_str(cvcp->wrap_center));\
	} \
	} \
	if ( ! cvcs.back()->name.size()) \
	cvcs.back()->name = std::string (def_config_key)+ \
	(cvm::to_str (++def_count)); \
	if (cvm::debug()) \
	cvm::log ("Done initializing a \""+ \
	std::string (def_config_key)+ \
	"\" component"+ \
	(cvm::debug() ? \
	", named \""+cvcs.back()->name+"\"" \
	: "")+".\n"); \
	def_conf = ""; \
	if (cvm::debug()) \
	cvm::log ("Parsed "+cvm::to_str (cvcs.size())+ \
	" components at this time.\n"); \
	} \
	}


	initialize_components ("distance", "distance", distance);
	initialize_components ("distance vector", "distanceVec", distance_vec);
	initialize_components ("distance vector "
	"direction", "distanceDir", distance_dir);
	initialize_components ("distance projection "
	"on an axis", "distanceZ", distance_z);
	initialize_components ("distance projection "
	"on a plane", "distanceXY", distance_xy);
	initialize_components ("minimum distance", "minDistance", min_distance);

	initialize_components ("coordination "
	"number", "coordNum", coordnum);
	initialize_components ("self-coordination "
	"number", "selfCoordNum", selfcoordnum);

	initialize_components ("angle", "angle", angle);
	initialize_components ("dihedral", "dihedral", dihedral);

	initialize_components ("hydrogen bond", "hBond", h_bond);

	// initialize_components ("alpha helix", "alphaDihedrals", alpha_dihedrals);
	initialize_components ("alpha helix", "alpha", alpha_angles);

	initialize_components ("dihedral principal "
	"component", "dihedralPC", dihedPC);

	initialize_components ("orientation", "orientation", orientation);
	initialize_components ("orientation "
	"angle", "orientationAngle",orientation_angle);
	initialize_components ("tilt", "tilt", tilt);
	initialize_components ("spin angle", "spinAngle", spin_angle);

	initialize_components ("RMSD", "rmsd", rmsd);

	// initialize_components ("logarithm of MSD", "logmsd", logmsd);

	initialize_components ("radius of "
	"gyration", "gyration", gyration);
	initialize_components ("eigenvector", "eigenvector", eigenvector);

	if (!cvcs.size())
	cvm::fatal_error ("Error: no valid components were provided "
	"for this collective variable.\n");

	cvm::log ("All components initialized.\n");


	// this is set false if any of the components has an exponent
	// different from 1 in the polynomial
	b_linear = true;

	// these will be set to false if any of the cvcs has them false
	b_inverse_gradients = true;
	b_Jacobian_force = true;

	// Decide whether the colvar is periodic
	// Used to wrap extended DOF if extendedLagrangian is on
	if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1
	&& (cvcs[0])->sup_coeff == 1.0 ) {
	this->b_periodic = true;
	this->period = (cvcs[0])->period;
	// TODO write explicit wrap() function for colvars to allow for
	// sup_coeff different from 1
	// this->period = (cvcs[0])->period * (cvcs[0])->sup_coeff;
	} else {
	this->b_periodic = false;
	this->period = 0.0;
	}

	// check the available features of each cvc
	for (size_t i = 0; i < cvcs.size(); i++) {

	if ((cvcs[i])->sup_np != 1) {
	if (cvm::debug() && b_linear)
	cvm::log ("Warning: You are using a non-linear polynomial "
	"combination to define this collective variable, "
	"some biasing methods may be unavailable.\n");
	b_linear = false;

	if ((cvcs[i])->sup_np < 0) {
	cvm::log ("Warning: you chose a negative exponent in the combination; "
	"if you apply forces, the simulation may become unstable "
	"when the component \""+
	(cvcs[i])->function_type+"\" approaches zero.\n");
	}
	}

	if ((cvcs[i])->b_periodic && !b_periodic) {
	cvm::log ("Warning: although this component is periodic, the colvar will "
	"not be treated as periodic, either because the exponent is not "
	"1, or because multiple components are present. Make sure that "
	"you know what you are doing!");
	}

	if (! (cvcs[i])->b_inverse_gradients)
	b_inverse_gradients = false;

	if (! (cvcs[i])->b_Jacobian_derivative)
	b_Jacobian_force = false;

	for (size_t j = i; j < cvcs.size(); j++) {
	if ( (cvcs[i])->type() != (cvcs[j])->type() ) {
	cvm::fatal_error ("ERROR: you are definining this collective variable "
	"by using components of different types, \""+
	colvarvalue::type_desc[(cvcs[i])->type()]+
	"\" and \""+
	colvarvalue::type_desc[(cvcs[j])->type()]+
	"\". "
	"You must use the same type in order to "
	" sum them together.\n");
	}
	}
	}


	{
	colvarvalue::Type const value_type = (cvcs[0])->type();
	if (cvm::debug())
	cvm::log ("This collective variable is a "+
	colvarvalue::type_desc[value_type]+", corresponding to "+
	cvm::to_str (colvarvalue::dof_num[value_type])+
	" internal degrees of freedom.\n");
	x.type (value_type);
	x_reported.type (value_type);
	}

	get_keyval (conf, "width", width, 1.0);
	if (width <= 0.0)
	cvm::fatal_error ("Error: \"width\" must be positive.\n");

	lower_boundary.type (this->type());
	lower_wall.type (this->type());

	upper_boundary.type (this->type());
	upper_wall.type (this->type());

	if (this->type() == colvarvalue::type_scalar) {

	if (get_keyval (conf, "lowerBoundary", lower_boundary, lower_boundary)) {
	enable (task_lower_boundary);
	}

	get_keyval (conf, "lowerWallConstant", lower_wall_k, 0.0);
	if (lower_wall_k > 0.0) {
	get_keyval (conf, "lowerWall", lower_wall, lower_boundary);
	enable (task_lower_wall);
	}

	if (get_keyval (conf, "upperBoundary", upper_boundary, upper_boundary)) {
	enable (task_upper_boundary);
	}

	get_keyval (conf, "upperWallConstant", upper_wall_k, 0.0);
	if (upper_wall_k > 0.0) {
	get_keyval (conf, "upperWall", upper_wall, upper_boundary);
	enable (task_upper_wall);
	}
	}

	// consistency checks for boundaries and walls
	if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) {
	if (lower_boundary >= upper_boundary) {
	cvm::fatal_error ("Error: the upper boundary, "+
	cvm::to_str (upper_boundary)+
	", is not higher than the lower boundary, "+
	cvm::to_str (lower_boundary)+".\n");
	}
	}

	if (tasks[task_lower_wall] && tasks[task_upper_wall]) {
	if (lower_wall >= upper_wall) {
	cvm::fatal_error ("Error: the upper wall, "+
	cvm::to_str (upper_wall)+
	", is not higher than the lower wall, "+
	cvm::to_str (lower_wall)+".\n");
	}

	if (dist2 (lower_wall, upper_wall) < 1.0E-12) {
	cvm::log ("Lower wall and upper wall are equal "
	"in the periodic domain of the colvar: disabling walls.\n");
	disable (task_lower_wall);
	disable (task_upper_wall);
	}
	}

	get_keyval (conf, "expandBoundaries", expand_boundaries, false);
	if (expand_boundaries && periodic_boundaries()) {
	cvm::fatal_error ("Error: trying to expand boundaries that already "
	"cover a whole period of a periodic colvar.\n");
	}

	{
	bool b_extended_lagrangian;
	get_keyval (conf, "extendedLagrangian", b_extended_lagrangian, false);

	if (b_extended_lagrangian) {
	cvm::real temp, tolerance, period;

	cvm::log ("Enabling the extended Lagrangian term for colvar \""+
	this->name+"\".\n");

	enable (task_extended_lagrangian);

	xr.type (this->type());
	vr.type (this->type());
	fr.type (this->type());

	const bool found = get_keyval (conf, "extendedTemp", temp, cvm::temperature());
	if (temp <= 0.0) {
	if (found)
	cvm::fatal_error ("Error: \"extendedTemp\" must be positive.\n");
	else
	cvm::fatal_error ("Error: a positive temperature must be provided, either "
	"by enabling a thermostat, or through \"extendedTemp\".\n");
	}

	get_keyval (conf, "extendedFluctuation", tolerance, 0.2*width);
	if (tolerance <= 0.0)
	cvm::fatal_error ("Error: \"extendedFluctuation\" must be positive.\n");
	ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
	cvm::log ("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");

	get_keyval (conf, "extendedTimeConstant", period, 40.0 * cvm::dt());
	if (period <= 0.0)
	cvm::fatal_error ("Error: \"extendedTimeConstant\" must be positive.\n");
	ext_mass = (cvm::boltzmann() * temp * period * period)
	/ (4.0 * PI * PI * tolerance * tolerance);
	cvm::log ("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)");

	{
	bool b_output_energy;
	get_keyval (conf, "outputEnergy", b_output_energy, false);
	if (b_output_energy) {
	enable (task_output_energy);
	}
	}

	get_keyval (conf, "extendedLangevinDamping", ext_gamma, 0.0);
	if (ext_gamma < 0.0)
	cvm::fatal_error ("Error: \"extendedLangevinDamping\" may not be negative.\n");
	if (ext_gamma != 0.0) {
	enable (task_langevin);
	ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1
	ext_sigma = sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / cvm::dt());
	}
	}
	}

	{
	bool b_output_value;
	get_keyval (conf, "outputValue", b_output_value, true);
	if (b_output_value) {
	enable (task_output_value);
	}
	}

	{
	bool b_output_velocity;
	get_keyval (conf, "outputVelocity", b_output_velocity, false);
	if (b_output_velocity) {
	enable (task_output_velocity);
	}
	}

	{
	bool b_output_system_force;
	get_keyval (conf, "outputSystemForce", b_output_system_force, false);
	if (b_output_system_force) {
	enable (task_output_system_force);
	}
	}

	{
	bool b_output_applied_force;
	get_keyval (conf, "outputAppliedForce", b_output_applied_force, false);
	if (b_output_applied_force) {
	enable (task_output_applied_force);
	}
	}

	if (cvm::b_analysis)
	parse_analysis (conf);

	if (cvm::debug())
	cvm::log ("Done initializing collective variable \""+this->name+"\".\n");
	}



	void colvar::build_atom_list (void)
	{
	// If atomic gradients are requested, build full list of atom ids from all cvcs
	std::list<int> temp_id_list;

	for (size_t i = 0; i < cvcs.size(); i++) {
	for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
	for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
	temp_id_list.push_back (cvcs[i]->atom_groups[j]->at(k).id);
	}
	}
	}

	temp_id_list.sort();
	temp_id_list.unique();
	atom_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
	temp_id_list.clear();

	atomic_gradients.resize (atom_ids.size());
	if (atom_ids.size()) {
	if (cvm::debug())
	cvm::log ("Colvar: created atom list with " + cvm::to_str(atom_ids.size()) + " atoms.\n");
	} else {
	cvm::log ("Warning: colvar components communicated no atom IDs.\n");
	}
	}


	void colvar::parse_analysis (std::string const &conf)
	{

	// if (cvm::debug())
	// cvm::log ("Parsing analysis flags for collective variable \""+
	// this->name+"\".\n");

	runave_length = 0;
	bool b_runave = false;
	if (get_keyval (conf, "runAve", b_runave) && b_runave) {

	enable (task_runave);

	get_keyval (conf, "runAveLength", runave_length, 1000);
	get_keyval (conf, "runAveStride", runave_stride, 1);

	if ((cvm::restart_out_freq % runave_stride) != 0)
	cvm::fatal_error ("Error: runAveStride must be commensurate with the restart frequency.\n");

	std::string runave_outfile;
	get_keyval (conf, "runAveOutputFile", runave_outfile,
	std::string (cvm::output_prefix+"."+
	this->name+".runave.traj"));

	size_t const this_cv_width = x.output_width (cvm::cv_width);
	runave_os.open (runave_outfile.c_str());
	runave_os << "# " << cvm::wrap_string ("step", cvm::it_width-2)
	<< " "
	<< cvm::wrap_string ("running average", this_cv_width)
	<< " "
	<< cvm::wrap_string ("running stddev", this_cv_width)
	<< "\n";
	}

	acf_length = 0;
	bool b_acf = false;
	if (get_keyval (conf, "corrFunc", b_acf) && b_acf) {

	enable (task_corrfunc);

	std::string acf_colvar_name;
	get_keyval (conf, "corrFuncWithColvar", acf_colvar_name, this->name);
	if (acf_colvar_name == this->name) {
	cvm::log ("Calculating auto-correlation function.\n");
	} else {
	cvm::log ("Calculating correlation function with \""+
	this->name+"\".\n");
	}

	std::string acf_type_str;
	get_keyval (conf, "corrFuncType", acf_type_str, to_lower_cppstr (std::string ("velocity")));
	if (acf_type_str == to_lower_cppstr (std::string ("coordinate"))) {
	acf_type = acf_coor;
	} else if (acf_type_str == to_lower_cppstr (std::string ("velocity"))) {
	acf_type = acf_vel;
	enable (task_fdiff_velocity);
	if (acf_colvar_name.size())
	(cvm::colvar_p (acf_colvar_name))->enable (task_fdiff_velocity);
	} else if (acf_type_str == to_lower_cppstr (std::string ("coordinate_p2"))) {
	acf_type = acf_p2coor;
	} else {
	cvm::fatal_error ("Unknown type of correlation function, \""+
	acf_type_str+"\".\n");
	}

	get_keyval (conf, "corrFuncOffset", acf_offset, 0);
	get_keyval (conf, "corrFuncLength", acf_length, 1000);
	get_keyval (conf, "corrFuncStride", acf_stride, 1);

	if ((cvm::restart_out_freq % acf_stride) != 0)
	cvm::fatal_error ("Error: corrFuncStride must be commensurate with the restart frequency.\n");

	get_keyval (conf, "corrFuncNormalize", acf_normalize, true);
	get_keyval (conf, "corrFuncOutputFile", acf_outfile,
	std::string (cvm::output_prefix+"."+this->name+
	".corrfunc.dat"));
	}
	}


	void colvar::enable (colvar::task const &t)
	{
	switch (t) {

	case task_output_system_force:
	enable (task_system_force);
	break;

	case task_report_Jacobian_force:
	enable (task_Jacobian_force);
	enable (task_system_force);
	if (cvm::debug())
	cvm::log ("Adding the Jacobian force to the system force, "
	"rather than applying its opposite silently.\n");
	break;

	case task_Jacobian_force:
	// checks below do not apply to extended-system colvars
	if ( !tasks[task_extended_lagrangian] ) {
	enable (task_gradients);

	if (!b_Jacobian_force)
	cvm::fatal_error ("Error: colvar \""+this->name+
	"\" does not have Jacobian forces implemented.\n");
	if (!b_linear)
	cvm::fatal_error ("Error: colvar \""+this->name+
	"\" must be defined as a linear combination "
	"to calculate the Jacobian force.\n");
	if (cvm::debug())
	cvm::log ("Enabling calculation of the Jacobian force "
	"on this colvar.\n");
	}
	fj.type (this->type());
	break;

	case task_system_force:
	if (!tasks[task_extended_lagrangian]) {
	if (!b_inverse_gradients) {
	cvm::fatal_error ("Error: one or more of the components of "
	"colvar \""+this->name+
	"\" does not support system force calculation.\n");
	}
	cvm::request_system_force();
	}
	ft.type (this->type());
	ft_reported.type (this->type());
	break;

	case task_output_applied_force:
	case task_lower_wall:
	case task_upper_wall:
	// all of the above require gradients
	enable (task_gradients);
	break;

	case task_fdiff_velocity:
	x_old.type (this->type());
	v_fdiff.type (this->type());
	v_reported.type (this->type());
	break;

	case task_output_velocity:
	enable (task_fdiff_velocity);
	break;

	case task_grid:
	if (this->type() != colvarvalue::type_scalar)
	cvm::fatal_error ("Cannot calculate a grid for collective variable, \""+
	this->name+"\", because its value is not a scalar number.\n");
	break;

	case task_extended_lagrangian:
	enable (task_gradients);
	v_reported.type (this->type());
	break;

	case task_lower_boundary:
	case task_upper_boundary:
	if (this->type() != colvarvalue::type_scalar) {
	cvm::fatal_error ("Error: this colvar is not a scalar value "
	"and cannot produce a grid.\n");
	}
	break;


	case task_output_value:
	case task_runave:
	case task_corrfunc:
	case task_ntot:
	break;

	case task_gradients:
	f.type (this->type());
	fb.type (this->type());
	break;

	case task_collect_gradients:
	if (this->type() != colvarvalue::type_scalar)
	cvm::fatal_error ("Collecting atomic gradients for non-scalar collective variable \""+
	this->name+"\" is not yet implemented.\n");
	enable (task_gradients);
	if (atom_ids.size() == 0) {
	build_atom_list();
	}
	break;
	}


	tasks[t] = true;
	}


	void colvar::disable (colvar::task const &t)
	{
	// check dependencies
	switch (t) {
	case task_gradients:
	disable (task_upper_wall);
	disable (task_lower_wall);
	disable (task_output_applied_force);
	disable (task_system_force);
	disable (task_Jacobian_force);
	break;

	case task_system_force:
	disable (task_output_system_force);
	break;

	case task_Jacobian_force:
	disable (task_report_Jacobian_force);
	break;

	case task_fdiff_velocity:
	disable (task_output_velocity);
	break;

	case task_lower_boundary:
	case task_upper_boundary:
	disable (task_grid);
	break;

	case task_extended_lagrangian:
	case task_report_Jacobian_force:
	case task_output_value:
	case task_output_velocity:
	case task_output_applied_force:
	case task_output_system_force:
	case task_runave:
	case task_corrfunc:
	case task_grid:
	case task_lower_wall:
	case task_upper_wall:
	case task_ntot:
	break;
	}

	tasks[t] = false;
	}


	colvar::~colvar()
	{
	if (cvm::b_analysis) {

	if (acf.size()) {
	cvm::log ("Writing acf to file \""+acf_outfile+"\".\n");

	std::ofstream acf_os (acf_outfile.c_str());
	if (! acf_os.good())
	cvm::fatal_error ("Cannot open file \""+acf_outfile+"\".\n");
	write_acf (acf_os);
	acf_os.close();
	}

	if (runave_os.good()) {
	runave_os.close();
	}
	}

	for (size_t i = 0; i < cvcs.size(); i++) {
	delete cvcs[i];
	}
	}



	// ****************** CALC FUNCTIONS ******************


	void colvar::calc()
	{
	if (cvm::debug())
	cvm::log ("Calculating colvar \""+this->name+"\".\n");

	// calculate the value of the colvar

	x.reset();
	if (x.type() == colvarvalue::type_scalar) {

	// scalar variable, polynomial combination allowed
	for (size_t i = 0; i < cvcs.size(); i++) {
	cvm::increase_depth();
	(cvcs[i])->calc_value();
	cvm::decrease_depth();
	if (cvm::debug())
	cvm::log ("Colvar component no. "+cvm::to_str (i+1)+
	" within colvar \""+this->name+"\" has value "+
	cvm::to_str ((cvcs[i])->value(),
	cvm::cv_width, cvm::cv_prec)+".\n");
	x += (cvcs[i])->sup_coeff *
	( ((cvcs[i])->sup_np != 1) ?
	std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
	(cvcs[i])->value().real_value );
	}
	} else {

	for (size_t i = 0; i < cvcs.size(); i++) {
	cvm::increase_depth();
	(cvcs[i])->calc_value();
	cvm::decrease_depth();
	if (cvm::debug())
	cvm::log ("Colvar component no. "+cvm::to_str (i+1)+
	" within colvar \""+this->name+"\" has value "+
	cvm::to_str ((cvcs[i])->value(),
	cvm::cv_width, cvm::cv_prec)+".\n");
	x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
	}
	}

	if (cvm::debug())
	cvm::log ("Colvar \""+this->name+"\" has value "+
	cvm::to_str (x, cvm::cv_width, cvm::cv_prec)+".\n");

	if (tasks[task_gradients]) {
	// calculate the gradients
	for (size_t i = 0; i < cvcs.size(); i++) {
	cvm::increase_depth();
	(cvcs[i])->calc_gradients();
	cvm::decrease_depth();
	}
	if (tasks[task_collect_gradients]) {
	// Collect the atomic gradients inside colvar object
	for (int a = 0; a < atomic_gradients.size(); a++) {
	atomic_gradients[a].reset();
	}
	for (size_t i = 0; i < cvcs.size(); i++) {
	// Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
	cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) *
	std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);

	for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {

	// If necessary, apply inverse rotation to get atomic
	// gradient in the laboratory frame
	if (cvcs[i]->atom_groups[j]->b_rotate) {
	cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse();

	for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
	int a = std::lower_bound (atom_ids.begin(), atom_ids.end(),
	cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
	atomic_gradients[a] += coeff *
	rot_inv.rotate (cvcs[i]->atom_groups[j]->at(k).grad);
	}

	} else {

	for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
	int a = std::lower_bound (atom_ids.begin(), atom_ids.end(),
	cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
	atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad;
	}
	}
	}
	}
	}
	}

	if (tasks[task_system_force]) {
	ft.reset();

	if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
	// get from the cvcs the system forces from the PREVIOUS step
	for (size_t i = 0; i < cvcs.size(); i++) {
	(cvcs[i])->calc_force_invgrads();
	// linear combination is assumed
	cvm::increase_depth();
	ft += (cvcs[i])->system_force() / ((cvcs[i])->sup_coeff * cvm::real (cvcs.size()));
	cvm::decrease_depth();
	}
	}

	if (tasks[task_report_Jacobian_force]) {
	// add the Jacobian force to the system force, and don't apply any silent
	// correction internally: biases such as colvarbias_abf will handle it
	ft += fj;
	}
	}

	if (tasks[task_fdiff_velocity]) {
	// calculate the velocity by finite differences
	if (cvm::step_relative() == 0)
	x_old = x;
	else {
	v_fdiff = fdiff_velocity (x_old, x);
	v_reported = v_fdiff;
	}
	}

	if (tasks[task_extended_lagrangian]) {

	// initialize the restraint center in the first step to the value
	// just calculated from the cvcs
	// TODO: put it in the restart information
	if (cvm::step_relative() == 0) {
	xr = x;
	vr = 0.0; // (already 0; added for clarity)
	}

	// report the restraint center as "value"
	x_reported = xr;
	v_reported = vr;
	// the "system force" with the extended Lagrangian is just the
	// harmonic term acting on the extended coordinate
	// Note: this is the force for current timestep
	ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad (xr, x);

	} else {

	x_reported = x;
	ft_reported = ft;
	}

	if (cvm::debug())
	cvm::log ("Done calculating colvar \""+this->name+"\".\n");
	}


	cvm::real colvar::update()
	{
	if (cvm::debug())
	cvm::log ("Updating colvar \""+this->name+"\".\n");


	// set to zero the applied force
	f.reset();

	// add the biases' force, which at this point should already have
	// been summed over each bias using this colvar
	f += fb;


	if (tasks[task_lower_wall] \|\| tasks[task_upper_wall]) {

	// wall force
	colvarvalue fw (this->type());

	// if the two walls are applied concurrently, decide which is the
	// closer one (on a periodic colvar, both walls may be applicable
	// at the same time)
	if ( (!tasks[task_upper_wall]) \|\|
	(this->dist2 (x, lower_wall) < this->dist2 (x, upper_wall)) ) {

	cvm::real const grad = this->dist2_lgrad (x, lower_wall);
	if (grad < 0.0) {
	fw = -0.5 * lower_wall_k * grad;
	if (cvm::debug())
	cvm::log ("Applying a lower wall force ("+
	cvm::to_str (fw)+") to \""+this->name+"\".\n");
	f += fw;

	}

	} else {

	cvm::real const grad = this->dist2_lgrad (x, upper_wall);
	if (grad > 0.0) {
	fw = -0.5 * upper_wall_k * grad;
	if (cvm::debug())
	cvm::log ("Applying an upper wall force ("+
	cvm::to_str (fw)+") to \""+this->name+"\".\n");
	f += fw;
	}
	}
	}

	if (tasks[task_Jacobian_force]) {

	cvm::increase_depth();
	for (size_t i = 0; i < cvcs.size(); i++) {
	(cvcs[i])->calc_Jacobian_derivative();
	}
	cvm::decrease_depth();

	fj.reset();
	for (size_t i = 0; i < cvcs.size(); i++) {
	// linear combination is assumed
	fj += 1.0 / ( cvm::real (cvcs.size()) * cvm::real ((cvcs[i])->sup_coeff) ) *
	(cvcs[i])->Jacobian_derivative();
	}
	fj = cvm::boltzmann() cvm::temperature();

	// the instantaneous Jacobian force was not included in the reported system force;
	// instead, it is subtracted from the applied force (silent Jacobian correction)
	if (! tasks[task_report_Jacobian_force])
	f -= fj;
	}


	if (tasks[task_extended_lagrangian]) {

	cvm::real dt = cvm::dt();

	// the total force is applied to the fictitious mass, while the
	// atoms only feel the harmonic force
	// fr: extended coordinate force; f: colvar force applied to atomic coordinates
	fr = f;
	fr += (-0.5 * ext_force_k) * this->dist2_lgrad (xr, x);
	f = (-0.5 * ext_force_k) * this->dist2_rgrad (xr, x);

	// leap frog: starting from x_i, f_i, v_(i-1/2)
	vr += (0.5 * dt) * fr / ext_mass;
	// Because of leapfrog, kinetic energy at time i is approximate
	kinetic_energy = 0.5 * ext_mass * vr * vr;
	potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
	// leap to v_(i+1/2)
	if (tasks[task_langevin]) {
	vr -= dt * ext_gamma * vr.real_value;
	vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
	}
	vr += (0.5 * dt) * fr / ext_mass;
	xr += dt * vr;
	xr.apply_constraints();
	if (this->b_periodic) this->wrap (xr);
	}


	if (tasks[task_fdiff_velocity]) {
	// set it for the next step
	x_old = x;
	}

	if (cvm::debug())
	cvm::log ("Done updating colvar \""+this->name+"\".\n");
	return (potential_energy + kinetic_energy);
	}


	void colvar::communicate_forces()
	{
	if (cvm::debug())
	cvm::log ("Communicating forces from colvar \""+this->name+"\".\n");

	if (x.type() == colvarvalue::type_scalar) {

	for (size_t i = 0; i < cvcs.size(); i++) {
	cvm::increase_depth();
	(cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff *
	cvm::real ((cvcs[i])->sup_np) *
	(std::pow ((cvcs[i])->value().real_value,
	(cvcs[i])->sup_np-1)) );
	cvm::decrease_depth();
	}

	} else {

	for (size_t i = 0; i < cvcs.size(); i++) {
	cvm::increase_depth();
	(cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff);
	cvm::decrease_depth();
	}
	}

	if (cvm::debug())
	cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n");
	}



	// ****************** METRIC FUNCTIONS ******************
	// Use the metrics defined by \link cvc \endlink objects


	bool colvar::periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) const
	{
	if ( (!tasks[task_lower_boundary]) \|\| (!tasks[task_upper_boundary]) ) {
	cvm::fatal_error ("Error: requesting to histogram the "
	"collective variable \""+this->name+"\", but a "
	"pair of lower and upper boundaries must be "
	"defined.\n");
	}

	if (period > 0.0) {
	if ( ((std::sqrt (this->dist2 (lb, ub))) / this->width)
	< 1.0E-10 ) {
	return true;
	}
	}

	return false;
	}

	bool colvar::periodic_boundaries() const
	{
	if ( (!tasks[task_lower_boundary]) \|\| (!tasks[task_upper_boundary]) ) {
	cvm::fatal_error ("Error: requesting to histogram the "
	"collective variable \""+this->name+"\", but a "
	"pair of lower and upper boundaries must be "
	"defined.\n");
	}

	return periodic_boundaries (lower_boundary, upper_boundary);
	}


	cvm::real colvar::dist2 (colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (cvcs[0])->dist2 (x1, x2);
	}

	colvarvalue colvar::dist2_lgrad (colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (cvcs[0])->dist2_lgrad (x1, x2);
	}

	colvarvalue colvar::dist2_rgrad (colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (cvcs[0])->dist2_rgrad (x1, x2);
	}

	cvm::real colvar::compare (colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (cvcs[0])->compare (x1, x2);
	}

	void colvar::wrap (colvarvalue &x) const
	{
	(cvcs[0])->wrap (x);
	return;
	}


	// ****************** INPUT FUNCTIONS ******************

	std::istream & colvar::read_restart (std::istream &is)
	{
	size_t const start_pos = is.tellg();

	std::string conf;
	if ( !(is >> colvarparse::read_block ("colvar", conf)) ) {
	// this is not a colvar block
	is.clear();
	is.seekg (start_pos, std::ios::beg);
	is.setstate (std::ios::failbit);
	return is;
	}

	{
	std::string check_name = "";
	if ( (get_keyval (conf, "name", check_name,
	std::string (""), colvarparse::parse_silent)) &&
	(check_name != name) ) {
	cvm::fatal_error ("Error: the state file does not match the "
	"configuration file, at colvar \""+name+"\".\n");
	}
	if (check_name.size() == 0) {
	cvm::fatal_error ("Error: Collective variable in the "
	"restart file without any identifier.\n");
	}
	}

	if ( !(get_keyval (conf, "x", x,
	colvarvalue (x.type()), colvarparse::parse_silent)) ) {
	cvm::log ("Error: restart file does not contain "
	"the value of the colvar \""+
	name+"\" .\n");
	} else {
	cvm::log ("Restarting collective variable \""+name+"\" from value: "+
	cvm::to_str (x)+"\n");
	}

	if (tasks[colvar::task_extended_lagrangian]) {

	if ( !(get_keyval (conf, "extended_x", xr,
	colvarvalue (x.type()), colvarparse::parse_silent)) &&
	!(get_keyval (conf, "extended_v", vr,
	colvarvalue (x.type()), colvarparse::parse_silent)) ) {
	cvm::log ("Error: restart file does not contain "
	"\"extended_x\" or \"extended_v\" for the colvar \""+
	name+"\", but you requested \"extendedLagrangian\".\n");
	}
	}

	if (tasks[task_extended_lagrangian]) {
	x_reported = xr;
	} else {
	x_reported = x;
	}

	if (tasks[task_output_velocity]) {

	if ( !(get_keyval (conf, "v", v_fdiff,
	colvarvalue (x.type()), colvarparse::parse_silent)) ) {
	cvm::log ("Error: restart file does not contain "
	"the velocity for the colvar \""+
	name+"\", but you requested \"outputVelocity\".\n");
	}

	if (tasks[task_extended_lagrangian]) {
	v_reported = vr;
	} else {
	v_reported = v_fdiff;
	}
	}

	return is;
	}


	std::istream & colvar::read_traj (std::istream &is)
	{
	size_t const start_pos = is.tellg();

	if (tasks[task_output_value]) {

	if (!(is >> x)) {
	cvm::log ("Error: in reading the value of colvar \""+
	this->name+"\" from trajectory.\n");
	is.clear();
	is.seekg (start_pos, std::ios::beg);
	is.setstate (std::ios::failbit);
	return is;
	}

	if (tasks[task_extended_lagrangian]) {
	is >> xr;
	x_reported = xr;
	} else {
	x_reported = x;
	}
	}

	if (tasks[task_output_velocity]) {

	is >> v_fdiff;

	if (tasks[task_extended_lagrangian]) {
	is >> vr;
	v_reported = vr;
	} else {
	v_reported = v_fdiff;
	}
	}

	if (tasks[task_output_system_force]) {

	is >> ft;

	if (tasks[task_extended_lagrangian]) {
	is >> fr;
	ft_reported = fr;
	} else {
	ft_reported = ft;
	}
	}

	if (tasks[task_output_applied_force]) {
	is >> f;
	}

	return is;
	}


	// ****************** OUTPUT FUNCTIONS ******************

	std::ostream & colvar::write_restart (std::ostream &os) {

	os << "colvar {\n"
	<< " name " << name << "\n"
	<< " x "
	<< std::setprecision (cvm::cv_prec)
	<< std::setw (cvm::cv_width)
	<< x << "\n";

	if (tasks[task_output_velocity]) {
	os << " v "
	<< std::setprecision (cvm::cv_prec)
	<< std::setw (cvm::cv_width)
	<< v_reported << "\n";
	}

	if (tasks[task_extended_lagrangian]) {
	os << " extended_x "
	<< std::setprecision (cvm::cv_prec)
	<< std::setw (cvm::cv_width)
	<< xr << "\n"
	<< " extended_v "
	<< std::setprecision (cvm::cv_prec)
	<< std::setw (cvm::cv_width)
	<< vr << "\n";
	}

	os << "}\n\n";

	return os;
	}


	std::ostream & colvar::write_traj_label (std::ostream & os)
	{
	size_t const this_cv_width = x.output_width (cvm::cv_width);

	os << " ";

	if (tasks[task_output_value]) {

	os << " "
	<< cvm::wrap_string (this->name, this_cv_width);

	if (tasks[task_extended_lagrangian]) {
	// restraint center
	os << " r_"
	<< cvm::wrap_string (this->name, this_cv_width-2);
	}
	}

	if (tasks[task_output_velocity]) {

	os << " v_"
	<< cvm::wrap_string (this->name, this_cv_width-2);

	if (tasks[task_extended_lagrangian]) {
	// restraint center
	os << " vr_"
	<< cvm::wrap_string (this->name, this_cv_width-3);
	}
	}

	if (tasks[task_output_energy]) {
	os << " Ep_"
	<< cvm::wrap_string (this->name, this_cv_width-3)
	<< " Ek_"
	<< cvm::wrap_string (this->name, this_cv_width-3);
	}

	if (tasks[task_output_system_force]) {

	os << " fs_"
	<< cvm::wrap_string (this->name, this_cv_width-2);

	if (tasks[task_extended_lagrangian]) {
	// restraint center
	os << " fr_"
	<< cvm::wrap_string (this->name, this_cv_width-3);
	}
	}

	if (tasks[task_output_applied_force]) {
	os << " fa_"
	<< cvm::wrap_string (this->name, this_cv_width-2);
	}

	return os;
	}


	std::ostream & colvar::write_traj (std::ostream &os)
	{
	os << " ";

	if (tasks[task_output_value]) {

	if (tasks[task_extended_lagrangian]) {
	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< x;
	}

	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< x_reported;
	}

	if (tasks[task_output_velocity]) {

	if (tasks[task_extended_lagrangian]) {
	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< v_fdiff;
	}

	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< v_reported;
	}

	if (tasks[task_output_energy]) {
	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< potential_energy
	<< " "
	<< kinetic_energy;
	}


	if (tasks[task_output_system_force]) {

	if (tasks[task_extended_lagrangian]) {
	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< ft;
	}

	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< ft_reported;
	}

	if (tasks[task_output_applied_force]) {
	os << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< f;
	}

	return os;
	}



	// ****************** ANALYSIS FUNCTIONS ******************

	void colvar::analyse()
	{
	if (tasks[task_runave]) {
	calc_runave();
	}

	if (tasks[task_corrfunc]) {
	calc_acf();
	}
	}


	inline void history_add_value (size_t const &history_length,
	std::list<colvarvalue> &history,
	colvarvalue const &new_value)
	{
	history.push_front (new_value);
	if (history.size() > history_length)
	history.pop_back();
	}

	inline void history_incr (std::list< std::list<colvarvalue> > &history,
	std::list< std::list<colvarvalue> >::iterator &history_p)
	{
	if ((++history_p) == history.end())
	history_p = history.begin();
	}


	void colvar::calc_acf()
	{
	// using here an acf_stride-long list of vectors for either
	// coordinates (acf_x_history) or velocities (acf_v_history); each vector can
	// contain up to acf_length values, which are contiguous in memory
	// representation but separated by acf_stride in the time series;
	// the pointer to each vector is changed at every step

	if (! (acf_x_history.size() \|\| acf_v_history.size()) ) {

	// first-step operations

	colvar *cfcv = (acf_colvar_name.size() ?
	cvm::colvar_p (acf_colvar_name) :
	this);
	if (cfcv->type() != this->type())
	cvm::fatal_error ("Error: correlation function between \""+cfcv->name+
	"\" and \""+this->name+"\" cannot be calculated, "
	"because their value types are different.\n");
	acf_nframes = 0;

	cvm::log ("Colvar \""+this->name+"\": initializing ACF calculation.\n");

	if (acf.size() < acf_length+1)
	acf.resize (acf_length+1, 0.0);

	switch (acf_type) {

	case acf_vel:
	// allocate space for the velocities history
	for (size_t i = 0; i < acf_stride; i++) {
	acf_v_history.push_back (std::list<colvarvalue>());
	}
	acf_v_history_p = acf_v_history.begin();
	break;

	case acf_coor:
	case acf_p2coor:
	// allocate space for the coordinates history
	for (size_t i = 0; i < acf_stride; i++) {
	acf_x_history.push_back (std::list<colvarvalue>());
	}
	acf_x_history_p = acf_x_history.begin();
	break;

	default:
	break;
	}

	} else {

	colvar *cfcv = (acf_colvar_name.size() ?
	cvm::colvar_p (acf_colvar_name) :
	this);

	switch (acf_type) {

	case acf_vel:

	if (tasks[task_fdiff_velocity]) {
	// calc() should do this already, but this only happens in a
	// simulation; better do it again in case a trajectory is
	// being read
	v_reported = v_fdiff = fdiff_velocity (x_old, cfcv->value());
	}

	calc_vel_acf ((*acf_v_history_p), cfcv->velocity());
	// store this value in the history
	history_add_value (acf_length+acf_offset, *acf_v_history_p, cfcv->velocity());
	// if stride is larger than one, cycle among different histories
	history_incr (acf_v_history, acf_v_history_p);
	break;

	case acf_coor:

	calc_coor_acf ((*acf_x_history_p), cfcv->value());
	history_add_value (acf_length+acf_offset, *acf_x_history_p, cfcv->value());
	history_incr (acf_x_history, acf_x_history_p);
	break;

	case acf_p2coor:

	calc_p2coor_acf ((*acf_x_history_p), cfcv->value());
	history_add_value (acf_length+acf_offset, *acf_x_history_p, cfcv->value());
	history_incr (acf_x_history, acf_x_history_p);
	break;

	default:
	break;
	}
	}

	if (tasks[task_fdiff_velocity]) {
	// set it for the next step
	x_old = x;
	}
	}


	void colvar::calc_vel_acf (std::list<colvarvalue> &v_list,
	colvarvalue const &v)
	{
	// loop over stored velocities and add to the ACF, but only the
	// length is sufficient to hold an entire row of ACF values
	if (v_list.size() >= acf_length+acf_offset) {
	std::list<colvarvalue>::iterator vs_i = v_list.begin();
	std::vector<cvm::real>::iterator acf_i = acf.begin();

	for (size_t i = 0; i < acf_offset; i++)
	vs_i++;

	// current vel with itself
	*(acf_i++) += v.norm2();

	// inner products of previous velocities with current (acf_i and
	// vs_i are updated)
	colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i);

	acf_nframes++;
	}
	}


	void colvar::calc_coor_acf (std::list<colvarvalue> &x_list,
	colvarvalue const &x)
	{
	// same as above but for coordinates
	if (x_list.size() >= acf_length+acf_offset) {
	std::list<colvarvalue>::iterator xs_i = x_list.begin();
	std::vector<cvm::real>::iterator acf_i = acf.begin();

	for (size_t i = 0; i < acf_offset; i++)
	xs_i++;

	*(acf_i++) += x.norm2();

	colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i);

	acf_nframes++;
	}
	}


	void colvar::calc_p2coor_acf (std::list<colvarvalue> &x_list,
	colvarvalue const &x)
	{
	// same as above but with second order Legendre polynomial instead
	// of just the scalar product
	if (x_list.size() >= acf_length+acf_offset) {
	std::list<colvarvalue>::iterator xs_i = x_list.begin();
	std::vector<cvm::real>::iterator acf_i = acf.begin();

	for (size_t i = 0; i < acf_offset; i++)
	xs_i++;

	// value of P2(0) = 1
	*(acf_i++) += 1.0;

	colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i);

	acf_nframes++;
	}
	}


	void colvar::write_acf (std::ostream &os)
	{
	if (!acf_nframes)
	cvm::log ("Warning: ACF was not calculated (insufficient frames).\n");
	os.setf (std::ios::scientific, std::ios::floatfield);
	os << "# Autocorrelation function for collective variable \""
	<< this->name << "\"\n";
	// one frame is used for normalization, the statistical sample is
	// hence decreased
	os << "# nframes = " << (acf_normalize ?
	acf_nframes - 1 :
	acf_nframes) << "\n";

	cvm::real const acf_norm = acf.front() / cvm::real (acf_nframes);
	std::vector<cvm::real>::iterator acf_i;
	size_t it = acf_offset;
	for (acf_i = acf.begin(); acf_i != acf.end(); acf_i++) {
	os << std::setw (cvm::it_width) << acf_stride * (it++) << " "
	<< std::setprecision (cvm::cv_prec)
	<< std::setw (cvm::cv_width)
	<< ( acf_normalize ?
	(acf_i)/(acf_norm cvm::real (acf_nframes)) :
	(*acf_i)/(cvm::real (acf_nframes)) ) << "\n";
	}
	}


	void colvar::calc_runave()
	{
	if (!x_history.size()) {

	runave.type (x.type());
	runave.reset();

	// first-step operations

	if (cvm::debug())
	cvm::log ("Colvar \""+this->name+
	"\": initializing running average calculation.\n");

	acf_nframes = 0;

	x_history.push_back (std::list<colvarvalue>());
	x_history_p = x_history.begin();

	} else {

	if ( (cvm::step_relative() % runave_stride) == 0) {

	if ((*x_history_p).size() >= runave_length-1) {

	runave = x;
	for (std::list<colvarvalue>::iterator xs_i = (*x_history_p).begin();
	xs_i != (*x_history_p).end(); xs_i++) {
	runave += (*xs_i);
	}
	runave *= 1.0 / cvm::real (runave_length);
	runave.apply_constraints();

	runave_variance = 0.0;
	runave_variance += this->dist2 (x, runave);
	for (std::list<colvarvalue>::iterator xs_i = (*x_history_p).begin();
	xs_i != (*x_history_p).end(); xs_i++) {
	runave_variance += this->dist2 (x, (*xs_i));
	}
	runave_variance *= 1.0 / cvm::real (runave_length-1);

	runave_os << std::setw (cvm::it_width) << cvm::step_relative()
	<< " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< runave << " "
	<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
	<< std::sqrt (runave_variance) << "\n";
	}

	history_add_value (runave_length, *x_history_p, x);
	}
	}

	}

colvar.CNo OneTemporaryActions

File Metadata

colvar.CView Options

Event Timeline

colvar.C
No OneTemporary
Actions

colvar.C
View Options