colvarcomp.h
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Mon, Nov 4, 00:04

colvarcomp.h
View Options

	// -- c++ --

	#ifndef COLVARCOMP_H
	#define COLVARCOMP_H

	// Declaration of colvar::cvc base class and derived ones.
	//
	// Future cvc's could be declared on additional header files.
	// After the declaration of a new derived class, its metric
	// functions must be reimplemented as well.
	// If the new cvc has no symmetry or periodicity,
	// this can be done straightforwardly by using the macro:
	// simple_scalar_dist_functions (derived_class)


	#include <fstream>
	#include <cmath>


	#include "colvarmodule.h"
	#include "colvar.h"
	#include "colvaratoms.h"


	/// \brief Colvar component (base class); most implementations of
	/// \link cvc \endlink utilize one or more \link
	/// colvarmodule::atom \endlink or \link colvarmodule::atom_group
	/// \endlink objects to access atoms.
	///
	/// A \link cvc \endlink object (or an object of a
	/// cvc-derived class) specifies how to calculate a collective
	/// variable, its gradients and other related physical quantities
	/// which do not depend only on the numeric value (the \link colvar
	/// \endlink class already serves this purpose).
	///
	/// No restriction is set to what kind of calculation a \link
	/// cvc \endlink object performs (usually calculate an
	/// analytical function of atomic coordinates). The only constraint
	/// is that the value calculated is implemented as a \link colvarvalue
	/// \endlink object. This serves to provide a unique way to calculate
	/// scalar and non-scalar collective variables, and specify if and how
	/// they can be combined together by the parent \link colvar \endlink
	/// object.
	///
	/// <b> If you wish to implement a new collective variable component, you
	/// should write your own class by inheriting directly from \link
	/// cvc \endlink, or one of its derived classes (for instance,
	/// \link distance \endlink is frequently used, because it provides
	/// useful data and function members for any colvar based on two
	/// atom groups). The steps are: \par
	/// 1. add the name of this class under colvar.h \par
	/// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par
	/// 3. declare the class in colvarcomp.h \par
	/// 4. implement the class in one of the files colvarcomp_*.C
	///
	/// </b>
	/// The cvm::atom and cvm::atom_group classes are available to
	/// transparently communicate with the simulation program. However,
	/// they are not strictly needed, as long as all the degrees of
	/// freedom associated to the cvc are properly evolved by a simple
	/// call to e.g. apply_force().

	class colvar::cvc
	: public colvarparse, public cvm::deps
	{
	public:

	/// \brief The name of the object (helps to identify this
	/// cvc instance when debugging)
	std::string name;

	/// \brief Description of the type of collective variable
	///
	/// Normally this string is set by the parent \link colvar \endlink
	/// object within its constructor, when all \link cvc \endlink
	/// objects are initialized; therefore the main "config string"
	/// constructor does not need to define it. If a \link cvc
	/// \endlink is initialized and/or a different constructor is used,
	/// this variable definition should be set within the constructor.
	std::string function_type;

	/// \brief Coefficient in the polynomial combination (default: 1.0)
	cvm::real sup_coeff;
	/// \brief Exponent in the polynomial combination (default: 1)
	int sup_np;

	/// \brief Is this a periodic component?
	bool b_periodic;

	/// \brief Period of this cvc value, (default: 0.0, non periodic)
	cvm::real period;

	/// \brief If the component is periodic, wrap around this value (default: 0.0)
	cvm::real wrap_center;

	/// \brief Constructor
	///
	/// At least one constructor which reads a string should be defined
	/// for every class inheriting from cvc \param conf Contents
	/// of the configuration file pertaining to this \link cvc
	/// \endlink
	cvc(std::string const &conf);

	/// \brief Within the constructor, make a group parse its own
	/// options from the provided configuration string
	/// Returns reference to new group
	cvm::atom_group *parse_group(std::string const &conf,
	char const *group_key,
	bool optional = false);

	/// \brief After construction, set data related to dependency handling
	int setup();

	/// \brief Default constructor (used when \link cvc \endlink
	/// objects are declared within other ones)
	cvc();

	/// Destructor
	virtual ~cvc();

	/// \brief Implementation of the feature list for colvar
	static std::vector<feature *> cvc_features;

	/// \brief Implementation of the feature list accessor for colvar
	virtual std::vector<feature *> &features() {
	return cvc_features;
	}

	/// \brief Obtain data needed for the calculation for the backend
	void read_data();

	/// \brief Calculate the variable
	virtual void calc_value() = 0;

	/// \brief Calculate the atomic gradients, to be reused later in
	/// order to apply forces
	virtual void calc_gradients() = 0;

	/// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
	virtual void debug_gradients(cvm::atom_group *group);

	/// \brief Calculate the total force from the system using the
	/// inverse atomic gradients
	virtual void calc_force_invgrads();

	/// \brief Calculate the divergence of the inverse atomic gradients
	virtual void calc_Jacobian_derivative();


	/// \brief Return the previously calculated value
	virtual colvarvalue const & value() const;

	// /// \brief Return const pointer to the previously calculated value
	// virtual const colvarvalue *p_value() const;

	/// \brief Return the previously calculated system force
	virtual colvarvalue const & system_force() const;

	/// \brief Return the previously calculated divergence of the
	/// inverse atomic gradients
	virtual colvarvalue const & Jacobian_derivative() const;

	/// \brief Apply the collective variable force, by communicating the
	/// atomic forces to the simulation program (\b Note: the \link ft
	/// \endlink member is not altered by this function)
	///
	/// Note: multiple calls to this function within the same simulation
	/// step will add the forces altogether \param cvforce The
	/// collective variable force, usually coming from the biases and
	/// eventually manipulated by the parent \link colvar \endlink
	/// object
	virtual void apply_force(colvarvalue const &cvforce) = 0;

	/// \brief Square distance between x1 and x2 (can be redefined to
	/// transparently implement constraints, symmetries and
	/// periodicities)
	///
	/// colvar::cvc::dist2() and the related functions are
	/// declared as "const" functions, but not "static", because
	/// additional parameters defining the metrics (e.g. the
	/// periodicity) may be specific to each colvar::cvc object.
	///
	/// If symmetries or periodicities are present, the
	/// colvar::cvc::dist2() should be redefined to return the
	/// "closest distance" value and colvar::cvc::dist2_lgrad(),
	/// colvar::cvc::dist2_rgrad() to return its gradients.
	///
	/// If constraints are present (and not already implemented by any
	/// of the \link colvarvalue \endlink types), the
	/// colvar::cvc::dist2_lgrad() and
	/// colvar::cvc::dist2_rgrad() functions should be redefined
	/// to provide a gradient which is compatible with the constraint,
	/// i.e. already deprived of its component normal to the constraint
	/// hypersurface.
	///
	/// Finally, another useful application, if you are performing very
	/// many operations with these functions, could be to override the
	/// \link colvarvalue \endlink member functions and access directly
	/// its member data. For instance: to define dist2(x1,x2) as
	/// (x2.real_value-x1.real_value)*(x2.real_value-x1.real_value) in
	/// case of a scalar \link colvarvalue \endlink type.
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;

	/// \brief Gradient(with respect to x1) of the square distance (can
	/// be redefined to transparently implement constraints, symmetries
	/// and periodicities)
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;

	/// \brief Gradient(with respect to x2) of the square distance (can
	/// be redefined to transparently implement constraints, symmetries
	/// and periodicities)
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;

	/// \brief Wrapp value (for periodic/symmetric cvcs)
	virtual void wrap(colvarvalue &x) const;

	/// \brief Pointers to all atom groups, to let colvars collect info
	/// e.g. atomic gradients
	std::vector<cvm::atom_group *> atom_groups;

	/// \brief Whether or not this CVC will be computed in parallel whenever possible
	bool b_try_scalable;

	protected:

	/// \brief Cached value
	colvarvalue x;

	/// \brief Value at the previous step
	colvarvalue x_old;

	/// \brief Calculated system force (\b Note: this is calculated from
	/// the total atomic forces read from the program, subtracting fromt
	/// the "internal" forces of the system the "external" forces from
	/// the colvar biases)
	colvarvalue ft;

	/// \brief Calculated Jacobian derivative (divergence of the inverse
	/// gradients): serves to calculate the phase space correction
	colvarvalue jd;
	};




	inline colvarvalue const & colvar::cvc::value() const
	{
	return x;
	}

	// inline const colvarvalue * colvar::cvc::p_value() const
	// {
	// return &x;
	// }

	inline colvarvalue const & colvar::cvc::system_force() const
	{
	return ft;
	}

	inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
	{
	return jd;
	}


	inline cvm::real colvar::cvc::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x1.dist2(x2);
	}

	inline colvarvalue colvar::cvc::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x1.dist2_grad(x2);
	}

	inline colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x2.dist2_grad(x1);
	}

	inline void colvar::cvc::wrap(colvarvalue &x) const
	{
	return;
	}


	/// \brief Colvar component: distance between the centers of mass of
	/// two groups (colvarvalue::type_scalar type, range [0:*))

	class colvar::distance
	: public colvar::cvc
	{
	protected:
	/// First atom group
	cvm::atom_group *group1;
	/// Second atom group
	cvm::atom_group *group2;
	/// Vector distance, cached to be recycled
	cvm::rvector dist_v;
	/// Use absolute positions, ignoring PBCs when present
	bool b_no_PBC;
	/// Compute system force on first site only to avoid unwanted
	/// coupling to other colvars (see e.g. Ciccotti et al., 2005)
	bool b_1site_force;
	public:
	distance(std::string const &conf);
	distance();
	virtual inline ~distance() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	// \brief Colvar component: distance vector between centers of mass
	// of two groups (\link colvarvalue::type_3vector \endlink type,
	// range (-:)x(-:)x(-:))
	class colvar::distance_vec
	: public colvar::distance
	{
	public:
	distance_vec(std::string const &conf);
	distance_vec();
	virtual inline ~distance_vec() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	/// Redefined to handle the box periodicity
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the box periodicity
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the box periodicity
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: distance unit vector (direction) between
	/// centers of mass of two groups (colvarvalue::type_unit3vector type,
	/// range [-1:1]x[-1:1]x[-1:1])
	class colvar::distance_dir
	: public colvar::distance
	{
	public:
	distance_dir(std::string const &conf);
	distance_dir();
	virtual inline ~distance_dir() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: projection of the distance vector along
	/// an axis(colvarvalue::type_scalar type, range (-:))
	class colvar::distance_z
	: public colvar::cvc
	{
	protected:
	/// Main atom group
	cvm::atom_group *main;
	/// Reference atom group
	cvm::atom_group *ref1;
	/// Optional, second ref atom group
	cvm::atom_group *ref2;
	/// Use absolute positions, ignoring PBCs when present
	bool b_no_PBC;
	/// Compute system force on one site only to avoid unwanted
	/// coupling to other colvars (see e.g. Ciccotti et al., 2005)
	bool b_1site_force;
	/// Vector on which the distance vector is projected
	cvm::rvector axis;
	/// Norm of the axis
	cvm::real axis_norm;
	/// Vector distance, cached to be recycled
	cvm::rvector dist_v;
	/// Flag: using a fixed axis vector?
	bool fixed_axis;
	public:
	distance_z(std::string const &conf);
	distance_z();
	virtual inline ~distance_z() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// \brief Redefined to make use of the user-provided period
	virtual void wrap(colvarvalue &x) const;
	};


	/// \brief Colvar component: projection of the distance vector on a
	/// plane (colvarvalue::type_scalar type, range [0:*))
	class colvar::distance_xy
	: public colvar::distance_z
	{
	protected:
	/// Components of the distance vector orthogonal to the axis
	cvm::rvector dist_v_ortho;
	/// Vector distances
	cvm::rvector v12, v13;
	public:
	distance_xy(std::string const &conf);
	distance_xy();
	virtual inline ~distance_xy() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
	/// as in NMR refinements(colvarvalue::type_scalar type, range (0:*))
	class colvar::distance_inv
	: public colvar::distance
	{
	protected:
	/// Components of the distance vector orthogonal to the axis
	int exponent;
	public:
	distance_inv(std::string const &conf);
	distance_inv();
	virtual inline ~distance_inv() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: N1xN2 vector of pairwise distances
	/// (colvarvalue::type_vector type, range (0:*) for each component)
	class colvar::distance_pairs
	: public colvar::cvc
	{
	protected:
	/// First atom group
	cvm::atom_group *group1;
	/// Second atom group
	cvm::atom_group *group2;
	/// Use absolute positions, ignoring PBCs when present
	bool b_no_PBC;
	public:
	distance_pairs(std::string const &conf);
	distance_pairs();
	virtual inline ~distance_pairs() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};



	/// \brief Colvar component: Radius of gyration of an atom group
	/// (colvarvalue::type_scalar type, range [0:*))
	class colvar::gyration
	: public colvar::cvc
	{
	protected:
	/// Atoms involved
	cvm::atom_group *atoms;
	public:
	/// Constructor
	gyration(std::string const &conf);
	gyration();
	virtual inline ~gyration() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: moment of inertia of an atom group
	/// (colvarvalue::type_scalar type, range [0:*))
	class colvar::inertia
	: public colvar::gyration
	{
	public:
	/// Constructor
	inertia(std::string const &conf);
	inertia();
	virtual inline ~inertia() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: moment of inertia of an atom group
	/// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
	class colvar::inertia_z
	: public colvar::inertia
	{
	protected:
	/// Vector on which the inertia tensor is projected
	cvm::rvector axis;
	public:
	/// Constructor
	inertia_z(std::string const &conf);
	inertia_z();
	virtual inline ~inertia_z() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: projection of 3N coordinates onto an
	/// eigenvector(colvarvalue::type_scalar type, range (-:))
	class colvar::eigenvector
	: public colvar::cvc
	{
	protected:

	/// Atom group
	cvm::atom_group * atoms;

	/// Reference coordinates
	std::vector<cvm::atom_pos> ref_pos;

	/// Geometric center of the reference coordinates
	cvm::atom_pos ref_pos_center;

	/// Eigenvector (of a normal or essential mode): will always have zero center
	std::vector<cvm::rvector> eigenvec;

	/// Inverse square norm of the eigenvector
	cvm::real eigenvec_invnorm2;

	public:

	/// Constructor
	eigenvector(std::string const &conf);
	virtual inline ~eigenvector() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};



	/// \brief Colvar component: angle between the centers of mass of
	/// three groups (colvarvalue::type_scalar type, range [0:PI])
	class colvar::angle
	: public colvar::cvc
	{
	protected:

	/// Atom group
	cvm::atom_group *group1;
	/// Atom group
	cvm::atom_group *group2;
	/// Atom group
	cvm::atom_group *group3;

	/// Inter site vectors
	cvm::rvector r21, r23;
	/// Inter site vector norms
	cvm::real r21l, r23l;
	/// Derivatives wrt group centers of mass
	cvm::rvector dxdr1, dxdr3;

	/// Compute system force on first site only to avoid unwanted
	/// coupling to other colvars (see e.g. Ciccotti et al., 2005)
	/// (or to allow dummy atoms)
	bool b_1site_force;
	public:

	/// Initialize by parsing the configuration
	angle(std::string const &conf);
	/// \brief Initialize the three groups after three atoms
	angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
	angle();
	virtual inline ~angle() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: angle between the dipole of a molecule and an axis
	/// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
	class colvar::dipole_angle
	: public colvar::cvc
	{
	protected:

	/// Dipole atom group
	cvm::atom_group *group1;
	/// Atom group
	cvm::atom_group *group2;
	/// Atom group
	cvm::atom_group *group3;

	/// Inter site vectors
	cvm::rvector r21, r23;
	/// Inter site vector norms
	cvm::real r21l, r23l;
	/// Derivatives wrt group centers of mass
	cvm::rvector dxdr1, dxdr3;

	/// Compute system force on first site only to avoid unwanted
	/// coupling to other colvars (see e.g. Ciccotti et al., 2005)
	/// (or to allow dummy atoms)
	bool b_1site_force;
	public:

	/// Initialize by parsing the configuration
	dipole_angle (std::string const &conf);
	/// \brief Initialize the three groups after three atoms
	dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
	dipole_angle();
	virtual inline ~dipole_angle() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force (colvarvalue const &force);
	virtual cvm::real dist2 (colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: dihedral between the centers of mass of
	/// four groups (colvarvalue::type_scalar type, range [-PI:PI])
	class colvar::dihedral
	: public colvar::cvc
	{
	protected:

	/// Atom group
	cvm::atom_group *group1;
	/// Atom group
	cvm::atom_group *group2;
	/// Atom group
	cvm::atom_group *group3;
	/// Atom group
	cvm::atom_group *group4;
	/// Inter site vectors
	cvm::rvector r12, r23, r34;

	/// \brief Compute system force on first site only to avoid unwanted
	/// coupling to other colvars (see e.g. Ciccotti et al., 2005)
	bool b_1site_force;

	public:

	/// Initialize by parsing the configuration
	dihedral(std::string const &conf);
	/// \brief Initialize the four groups after four atoms
	dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
	dihedral();
	virtual inline ~dihedral() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);

	/// Redefined to handle the 2*PI periodicity
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual void wrap(colvarvalue &x) const;
	};


	/// \brief Colvar component: coordination number between two groups
	/// (colvarvalue::type_scalar type, range [0:N1*N2])
	class colvar::coordnum
	: public colvar::cvc
	{
	protected:
	/// First atom group
	cvm::atom_group *group1;
	/// Second atom group
	cvm::atom_group *group2;
	/// \brief "Cutoff" for isotropic calculation (default)
	cvm::real r0;
	/// \brief "Cutoff vector" for anisotropic calculation
	cvm::rvector r0_vec;
	/// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
	/// used
	bool b_anisotropic;
	/// Integer exponent of the function numerator
	int en;
	/// Integer exponent of the function denominator
	int ed;
	/// \brief If true, group2 will be treated as a single atom
	/// (default: loop over all pairs of atoms in group1 and group2)
	bool b_group2_center_only;
	public:
	/// Constructor
	coordnum(std::string const &conf);
	coordnum();
	virtual inline ~coordnum() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	template<bool b_gradients>
	/// \brief Calculate a coordination number through the function
	/// (1-xn)/(1-xm), x = \|A1-A2\|/r0 \param r0 "cutoff" for the
	/// coordination number \param exp_num \i n exponent \param exp_den
	/// \i m exponent \param A1 atom \param A2 atom
	static cvm::real switching_function(cvm::real const &r0,
	int const &exp_num, int const &exp_den,
	cvm::atom &A1, cvm::atom &A2);

	template<bool b_gradients>
	/// \brief Calculate a coordination number through the function
	/// (1-xn)/(1-xm), x = \|(A1-A2)*(r0_vec)^-\|1 \param r0_vec
	/// vector of different cutoffs in the three directions \param
	/// exp_num \i n exponent \param exp_den \i m exponent \param A1
	/// atom \param A2 atom
	static cvm::real switching_function(cvm::rvector const &r0_vec,
	int const &exp_num, int const &exp_den,
	cvm::atom &A1, cvm::atom &A2);

	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: self-coordination number within a group
	/// (colvarvalue::type_scalar type, range [0:N*(N-1)/2])
	class colvar::selfcoordnum
	: public colvar::cvc
	{
	protected:
	/// First atom group
	cvm::atom_group *group1;
	/// \brief "Cutoff" for isotropic calculation (default)
	cvm::real r0;
	/// Integer exponent of the function numerator
	int en;
	/// Integer exponent of the function denominator
	int ed;
	public:
	/// Constructor
	selfcoordnum(std::string const &conf);
	selfcoordnum();
	virtual inline ~selfcoordnum() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	template<bool b_gradients>
	/// \brief Calculate a coordination number through the function
	/// (1-xn)/(1-xm), x = \|A1-A2\|/r0 \param r0 "cutoff" for the
	/// coordination number \param exp_num \i n exponent \param exp_den
	/// \i m exponent \param A1 atom \param A2 atom
	static cvm::real switching_function(cvm::real const &r0,
	int const &exp_num, int const &exp_den,
	cvm::atom &A1, cvm::atom &A2);

	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: hydrogen bond, defined as the product of
	/// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
	/// (colvarvalue::type_scalar type, range [0:1])
	class colvar::h_bond
	: public colvar::cvc
	{
	protected:
	/// \brief "Cutoff" distance between acceptor and donor
	cvm::real r0;
	/// Integer exponent of the function numerator
	int en;
	/// Integer exponent of the function denominator
	int ed;
	public:
	h_bond(std::string const &conf);
	/// Constructor for atoms already allocated
	h_bond(cvm::atom const &acceptor,
	cvm::atom const &donor,
	cvm::real r0, int en, int ed);
	h_bond();
	virtual ~h_bond();
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);

	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: alpha helix content of a contiguous
	/// segment of 5 or more residues, implemented as a sum of phi/psi
	/// dihedral angles and hydrogen bonds (colvarvalue::type_scalar type,
	/// range [0:1])
	// class colvar::alpha_dihedrals
	// : public colvar::cvc
	// {
	// protected:

	// /// Alpha-helical reference phi value
	// cvm::real phi_ref;

	// /// Alpha-helical reference psi value
	// cvm::real psi_ref;

	// /// List of phi dihedral angles
	// std::vector<dihedral *> phi;

	// /// List of psi dihedral angles
	// std::vector<dihedral *> psi;

	// /// List of hydrogen bonds
	// std::vector<h_bond *> hb;

	// public:

	// alpha_dihedrals (std::string const &conf);
	// alpha_dihedrals();
	// virtual inline ~alpha_dihedrals() {}
	// virtual void calc_value();
	// virtual void calc_gradients();
	// virtual void apply_force (colvarvalue const &force);
	// virtual cvm::real dist2 (colvarvalue const &x1,
	// colvarvalue const &x2) const;
	// virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
	// colvarvalue const &x2) const;
	// virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
	// colvarvalue const &x2) const;
	// };


	/// \brief Colvar component: alpha helix content of a contiguous
	/// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca
	/// angles and hydrogen bonds (colvarvalue::type_scalar type, range
	/// [0:1])
	class colvar::alpha_angles
	: public colvar::cvc
	{
	protected:

	/// Reference Calpha-Calpha angle (default: 88 degrees)
	cvm::real theta_ref;

	/// Tolerance on the Calpha-Calpha angle
	cvm::real theta_tol;

	/// List of Calpha-Calpha angles
	std::vector<angle *> theta;

	/// List of hydrogen bonds
	std::vector<h_bond *> hb;

	/// Contribution of the hb terms
	cvm::real hb_coeff;

	public:

	alpha_angles(std::string const &conf);
	alpha_angles();
	virtual ~alpha_angles();
	void calc_value();
	void calc_gradients();
	void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: dihedPC
	/// Projection of the config onto a dihedral principal component
	/// See e.g. Altis et al., J. Chem. Phys 126, 244111 (2007)
	/// Based on a set of 'dihedral' cvcs
	class colvar::dihedPC
	: public colvar::cvc
	{
	protected:

	std::vector<dihedral *> theta;
	std::vector<cvm::real> coeffs;

	public:

	dihedPC(std::string const &conf);
	dihedPC();
	virtual ~dihedPC();
	void calc_value();
	void calc_gradients();
	void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};

	/// \brief Colvar component: orientation in space of an atom group,
	/// with respect to a set of reference coordinates
	/// (colvarvalue::type_quaternion type, range
	/// [-1:1]x[-1:1]x[-1:1]x[-1:1])
	class colvar::orientation
	: public colvar::cvc
	{
	protected:

	/// Atom group
	cvm::atom_group * atoms;
	/// Center of geometry of the group
	cvm::atom_pos atoms_cog;

	/// Reference coordinates
	std::vector<cvm::atom_pos> ref_pos;

	/// Rotation object
	cvm::rotation rot;

	/// \brief This is used to remove jumps in the sign of the
	/// quaternion, which may be annoying in the colvars trajectory
	cvm::quaternion ref_quat;

	public:

	orientation(std::string const &conf);
	orientation();
	virtual inline ~orientation() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: angle of rotation with respect to a set
	/// of reference coordinates (colvarvalue::type_scalar type, range
	/// [0:PI))
	class colvar::orientation_angle
	: public colvar::orientation
	{
	public:

	orientation_angle(std::string const &conf);
	orientation_angle();
	virtual inline ~orientation_angle() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: cosine of the angle of rotation with respect to a set
	/// of reference coordinates (colvarvalue::type_scalar type, range
	/// [-1:1])
	class colvar::orientation_proj
	: public colvar::orientation
	{
	public:

	orientation_proj(std::string const &conf);
	orientation_proj();
	virtual inline ~orientation_proj() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	/// \brief Colvar component: projection of the orientation vector onto
	/// a predefined axis (colvarvalue::type_scalar type, range [-1:1])
	class colvar::tilt
	: public colvar::orientation
	{
	protected:

	cvm::rvector axis;

	public:

	tilt(std::string const &conf);
	tilt();
	virtual inline ~tilt() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};



	/// \brief Colvar component: angle of rotation around a predefined
	/// axis (colvarvalue::type_scalar type, range [-PI:PI])
	class colvar::spin_angle
	: public colvar::orientation
	{
	protected:

	cvm::rvector axis;

	public:

	spin_angle(std::string const &conf);
	spin_angle();
	virtual inline ~spin_angle() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	/// Redefined to handle the 2*PI periodicity
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	/// Redefined to handle the 2*PI periodicity
	virtual void wrap(colvarvalue &x) const;
	};



	/// \brief Colvar component: root mean square deviation (RMSD) of a
	/// group with respect to a set of reference coordinates; uses \link
	/// colvar::orientation \endlink to calculate the rotation matrix
	/// (colvarvalue::type_scalar type, range [0:*))
	class colvar::rmsd
	: public colvar::cvc
	{
	protected:

	/// Atom group
	cvm::atom_group *atoms;

	/// Reference coordinates (for RMSD calculation only)
	std::vector<cvm::atom_pos> ref_pos;

	public:

	/// Constructor
	rmsd(std::string const &conf);
	virtual inline ~rmsd() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void calc_force_invgrads();
	virtual void calc_Jacobian_derivative();
	virtual void apply_force(colvarvalue const &force);
	virtual cvm::real dist2(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const;
	};


	// \brief Colvar component: flat vector of Cartesian coordinates
	// Mostly useful to compute scripted colvar values
	class colvar::cartesian
	: public colvar::cvc
	{
	protected:
	/// Atom group
	cvm::atom_group *atoms;
	/// Which Cartesian coordinates to include
	std::vector<size_t> axes;
	public:
	cartesian(std::string const &conf);
	cartesian();
	virtual inline ~cartesian() {}
	virtual void calc_value();
	virtual void calc_gradients();
	virtual void apply_force(colvarvalue const &force);
	};


	// metrics functions for cvc implementations

	// simple definitions of the distance functions; these are useful only
	// for optimization (the type check performed in the default
	// colvarcomp functions is skipped)

	// definitions assuming the scalar type

	#define simple_scalar_dist_functions(TYPE) \
	\
	inline cvm::real colvar::TYPE::dist2(colvarvalue const &x1, \
	colvarvalue const &x2) const \
	{ \
	return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
	} \
	\
	inline colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1, \
	colvarvalue const &x2) const \
	{ \
	return 2.0 * (x1.real_value - x2.real_value); \
	} \
	\
	inline colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1, \
	colvarvalue const &x2) const \
	{ \
	return this->dist2_lgrad(x2, x1); \
	} \
	\

	simple_scalar_dist_functions(distance)
	// NOTE: distance_z has explicit functions, see below
	simple_scalar_dist_functions(distance_xy)
	simple_scalar_dist_functions(distance_inv)
	simple_scalar_dist_functions(angle)
	simple_scalar_dist_functions(dipole_angle)
	simple_scalar_dist_functions(coordnum)
	simple_scalar_dist_functions(selfcoordnum)
	simple_scalar_dist_functions(h_bond)
	simple_scalar_dist_functions(gyration)
	simple_scalar_dist_functions(inertia)
	simple_scalar_dist_functions(inertia_z)
	simple_scalar_dist_functions(rmsd)
	simple_scalar_dist_functions(orientation_angle)
	simple_scalar_dist_functions(orientation_proj)
	simple_scalar_dist_functions(tilt)
	simple_scalar_dist_functions(eigenvector)
	// simple_scalar_dist_functions (alpha_dihedrals)
	simple_scalar_dist_functions(alpha_angles)
	simple_scalar_dist_functions(dihedPC)


	// metrics functions for cvc implementations with a periodicity

	inline cvm::real colvar::dihedral::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return diff * diff;
	}

	inline colvarvalue colvar::dihedral::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return 2.0 * diff;
	}

	inline colvarvalue colvar::dihedral::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return (-2.0) * diff;
	}

	inline void colvar::dihedral::wrap(colvarvalue &x) const
	{
	if ((x.real_value - wrap_center) >= 180.0) {
	x.real_value -= 360.0;
	return;
	}

	if ((x.real_value - wrap_center) < -180.0) {
	x.real_value += 360.0;
	return;
	}

	return;
	}

	inline cvm::real colvar::spin_angle::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return diff * diff;
	}

	inline colvarvalue colvar::spin_angle::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return 2.0 * diff;
	}

	inline colvarvalue colvar::spin_angle::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
	return (-2.0) * diff;
	}

	inline void colvar::spin_angle::wrap(colvarvalue &x) const
	{
	if ((x.real_value - wrap_center) >= 180.0) {
	x.real_value -= 360.0;
	return;
	}

	if ((x.real_value - wrap_center) < -180.0) {
	x.real_value += 360.0;
	return;
	}

	return;
	}


	// Projected distance
	// Differences should always be wrapped around 0 (ignoring wrap_center)
	inline cvm::real colvar::distance_z::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	if (period != 0.0) {
	cvm::real shift = std::floor(diff/period + 0.5);
	diff -= shift * period;
	}
	return diff * diff;
	}

	inline colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	if (period != 0.0) {
	cvm::real shift = std::floor(diff/period + 0.5);
	diff -= shift * period;
	}
	return 2.0 * diff;
	}

	inline colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	cvm::real diff = x1.real_value - x2.real_value;
	if (period != 0.0) {
	cvm::real shift = std::floor(diff/period + 0.5);
	diff -= shift * period;
	}
	return (-2.0) * diff;
	}

	inline void colvar::distance_z::wrap(colvarvalue &x) const
	{
	if (! this->b_periodic) {
	// don't wrap if the period has not been set
	return;
	}

	cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
	x.real_value -= shift * period;
	return;
	}


	// distance between three dimensional vectors
	//
	// TODO apply PBC to distance_vec
	// Note: differences should be centered around (0, 0, 0)!

	inline cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
	}

	inline colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
	}

	inline colvarvalue colvar::distance_vec::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
	}

	inline cvm::real colvar::distance_dir::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (x1.rvector_value - x2.rvector_value).norm2();
	}

	inline colvarvalue colvar::distance_dir::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return colvarvalue((x1.rvector_value - x2.rvector_value), colvarvalue::type_unit3vector);
	}

	inline colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector);
	}

	inline cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return (x1.vector1d_value - x2.vector1d_value).norm2();
	}

	inline colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return colvarvalue((x1.vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
	}

	inline colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return colvarvalue((x2.vector1d_value - x1.vector1d_value), colvarvalue::type_vector);
	}


	// distance between quaternions

	inline cvm::real colvar::orientation::dist2(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x1.quaternion_value.dist2(x2);
	}

	inline colvarvalue colvar::orientation::dist2_lgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x1.quaternion_value.dist2_grad(x2);
	}

	inline colvarvalue colvar::orientation::dist2_rgrad(colvarvalue const &x1,
	colvarvalue const &x2) const
	{
	return x2.quaternion_value.dist2_grad(x1);
	}


	#endif

colvarcomp.hNo OneTemporaryActions

File Metadata

colvarcomp.hView Options

Event Timeline

colvarcomp.h
No OneTemporary
Actions

colvarcomp.h
View Options