Diffusion lammps (efficient_neuronet) Browse

colvarcomp.h
da8b07b4f090efficient_neuronet

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

lammps/lib/colvars/colvarcomp.hda8b07b4f090efficient_neuronet

colvarcomp.hView Options

lammps/lib/colvars/colvarcomp.h
da8b07b4f090efficient_neuronet

colvarcomp.h
View Options