Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85648276
colvar.h
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Sep 30, 17:17
Size
18 KB
Mime Type
text/x-c++
Expires
Wed, Oct 2, 17:17 (2 d)
Engine
blob
Format
Raw Data
Handle
21227115
Attached To
rLAMMPS lammps
colvar.h
View Options
// -*- c++ -*-
#ifndef COLVAR_H
#define COLVAR_H
#include <iostream>
#include <iomanip>
#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvardeps.h"
/// \brief A collective variable (main class); to be defined, it needs
/// at least one object of a derived class of colvar::cvc; it
/// calculates and returns a \link colvarvalue \endlink object
///
/// This class parses the configuration, defines the behaviour and
/// stores the value (\link colvar::x \endlink) and all related data
/// of a collective variable. How the value is calculated is defined
/// in \link colvar::cvc \endlink and its derived classes. The
/// \link colvar \endlink object contains pointers to multiple \link
/// colvar::cvc \endlink derived objects, which can be combined
/// together into one collective variable. This makes possible to
/// implement new collective variables at runtime based on the
/// existing ones. Currently, this possibility is limited to a
/// polynomial, using the coefficients cvc::sup_coeff and the
/// exponents cvc::sup_np. In case of non-scalar variables,
/// only exponents equal to 1 are accepted.
///
/// Please note that most of its members are \link colvarvalue
/// \endlink objects, i.e. they can handle different data types
/// together, and must all be set to the same type of colvar::value()
/// before using them together in assignments or other operations; this is usually done
/// automatically in the constructor. If you add a new member of
/// \link colvarvalue \endlink type, you should also add its
/// initialization line in the \link colvar \endlink constructor.
class colvar : public colvarparse, public cvm::deps {
public:
/// Name
std::string name;
/// \brief Current value (previously set by calc() or by read_traj())
colvarvalue const & value() const;
/// \brief Current actual value (not extended DOF)
colvarvalue const & actual_value() const;
/// \brief Current velocity (previously set by calc() or by read_traj())
colvarvalue const & velocity() const;
/// \brief Current system force (previously obtained from calc() or
/// read_traj()). Note: this is calculated using the atomic forces
/// from the last simulation step.
///
/// Total atom forces are read from the MD program, the total force
/// acting on the collective variable is calculated summing those
/// from all colvar components, the bias and walls forces are
/// subtracted.
colvarvalue const & system_force() const;
/// \brief Typical fluctuation amplitude for this collective
/// variable (e.g. local width of a free energy basin)
///
/// In metadynamics calculations, \link colvarbias_meta \endlink,
/// this value is used to calculate the width of a gaussian. In ABF
/// calculations, \link colvarbias_abf \endlink, it is used to
/// calculate the grid spacing in the direction of this collective
/// variable.
cvm::real width;
/// \brief Implementation of the feature list for colvar
static std::vector<feature *> cv_features;
/// \brief Implementation of the feature list accessor for colvar
std::vector<feature *> &features() {
return cv_features;
}
int refresh_deps();
/// List of biases that depend on this colvar
std::vector<colvarbias *> biases;
protected:
/*
extended:
H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
+ \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
+ \sum_{t'<t} W * exp(-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
+ \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M
normal:
H = H_{0} + \sum_{t'<t} W * exp(-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
+ \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M
output:
H = H_{0} (only output S(x), no forces)
Here:
S(x(t)) = x
s(t) = xr
DS = Ds = delta
*/
/// Value of the colvar
colvarvalue x;
// TODO: implement functionality to treat these
// /// Vector of individual values from CVCs
// colvarvalue x_cvc;
// /// Jacobian matrix of individual values from CVCs
// colvarvalue dx_cvc;
/// Cached reported value (x may be manipulated)
colvarvalue x_reported;
/// Finite-difference velocity
colvarvalue v_fdiff;
inline colvarvalue fdiff_velocity(colvarvalue const &xold,
colvarvalue const &xnew)
{
// using the gradient of the square distance to calculate the
// velocity (non-scalar variables automatically taken into
// account)
cvm::real dt = cvm::dt();
return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) *
0.5 * dist2_lgrad(xnew, xold) );
}
/// Cached reported velocity
colvarvalue v_reported;
// Options for extended_lagrangian
/// Restraint center
colvarvalue xr;
/// Velocity of the restraint center
colvarvalue vr;
/// Mass of the restraint center
cvm::real ext_mass;
/// Restraint force constant
cvm::real ext_force_k;
/// Friction coefficient for Langevin extended dynamics
cvm::real ext_gamma;
/// Amplitude of Gaussian white noise for Langevin extended dynamics
cvm::real ext_sigma;
/// \brief Harmonic restraint force
colvarvalue fr;
/// \brief Jacobian force, when Jacobian_force is enabled
colvarvalue fj;
/// Cached reported system force
colvarvalue ft_reported;
public:
/// \brief Bias force; reset_bias_force() should be called before
/// the biases are updated
colvarvalue fb;
/// \brief Total \em applied force; fr (if extended_lagrangian
/// is defined), fb (if biases are applied) and the walls' forces
/// (if defined) contribute to it
colvarvalue f;
/// \brief Total force, as derived from the atomic trajectory;
/// should equal the total system force plus \link f \endlink
colvarvalue ft;
/// Period, if it is a constant
cvm::real period;
/// \brief Same as above, but also takes into account components
/// with a variable period, such as distanceZ
bool b_periodic;
/// \brief Expand the boundaries of multiples of width, to keep the
/// value always within range
bool expand_boundaries;
/// \brief Location of the lower boundary
colvarvalue lower_boundary;
/// \brief Location of the lower wall
colvarvalue lower_wall;
/// \brief Force constant for the lower boundary potential (|x-xb|^2)
cvm::real lower_wall_k;
/// \brief Whether this colvar has a hard lower boundary
bool hard_lower_boundary;
/// \brief Location of the upper boundary
colvarvalue upper_boundary;
/// \brief Location of the upper wall
colvarvalue upper_wall;
/// \brief Force constant for the upper boundary potential (|x-xb|^2)
cvm::real upper_wall_k;
/// \brief Whether this colvar has a hard upper boundary
bool hard_upper_boundary;
/// \brief Is the interval defined by the two boundaries periodic?
bool periodic_boundaries() const;
/// \brief Is the interval defined by the two boundaries periodic?
bool periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const;
/// Constructor
colvar(std::string const &conf);
/// Parse the CVC configuration and allocate their data
int init_components(std::string const &conf);
private:
/// Parse the CVC configuration for all components of a certain type
template<typename def_class_name> int init_components_type(std::string const &conf,
char const *def_desc,
char const *def_config_key);
public:
/// Get ready for a run and re-initialize internal data if needed
void setup();
/// Destructor
~colvar();
/// \brief Calculate the colvar's value and related quantities
int calc();
/// \brief Calculate a subset of the colvar components (CVCs) currently active
/// (default: all active CVCs)
/// Note: both arguments refer to the sect of *active* CVCs, not all CVCs
int calc_cvcs(int first = 0, size_t num_cvcs = 0);
/// Ensure that the selected range of CVCs is consistent
int check_cvc_range(int first_cvc, size_t num_cvcs);
/// \brief Calculate the values of the given subset of CVCs
int calc_cvc_values(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for gradients
int calc_cvc_gradients(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for system forces
int calc_cvc_sys_forces(int first, size_t num_cvcs);
/// \brief Same as \link colvar::calc_cvc_values \endlink but for Jacobian derivatives/forces
int calc_cvc_Jacobians(int first, size_t num_cvcs);
/// \brief Collect quantities from CVCs and update aggregated data for the colvar
int collect_cvc_data();
/// \brief Collect the values of the CVCs
int collect_cvc_values();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for gradients
int collect_cvc_gradients();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for system forces
int collect_cvc_sys_forces();
/// \brief Same as \link colvar::collect_cvc_values \endlink but for Jacobian derivatives/forces
int collect_cvc_Jacobians();
/// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
int calc_colvar_properties();
/// Get the current biasing force
inline colvarvalue bias_force() const
{
return fb;
}
/// Set the total biasing force to zero
void reset_bias_force();
/// Add to the total force from biases
void add_bias_force(colvarvalue const &force);
/// \brief Collect all forces on this colvar, integrate internal
/// equations of motion of internal degrees of freedom; see also
/// colvar::communicate_forces()
/// return colvar energy if extended Lagrandian active
cvm::real update_forces_energy();
/// \brief Communicate forces (previously calculated in
/// colvar::update()) to the external degrees of freedom
void communicate_forces();
/// \brief Enables and disables individual CVCs based on flags
int set_cvc_flags(std::vector<bool> const &flags);
/// \brief Updates the flags in the CVC objects, and their number
int update_cvc_flags();
protected:
/// \brief Number of CVC objects with an active flag
size_t n_active_cvcs;
/// Sum of square coefficients for active cvcs
cvm::real active_cvc_square_norm;
/// Time step multiplier (for coarse-time-step colvars)
/// Colvar will only be calculated at those times; biases may ignore the information and
/// always update their own forces (which is typically inexpensive) especially if
/// they rely on other colvars. In this case, the colvar will accumulate forces applied between
/// colvar updates. Alternately they may use it to calculate "impulse" biasing
/// forces at longer intervals. Impulse forces must be multiplied by the timestep factor.
int time_step_factor;
/// Biasing force collected between updates, to be applied at next update for coarse-time-step colvars
colvarvalue f_accumulated;
public:
/// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
inline size_t num_active_cvcs() const { return n_active_cvcs; }
/// \brief returns time_step_factor
inline int get_time_step_factor() const {return time_step_factor;}
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval
///
/// Handles correctly symmetries and periodic boundary conditions
void wrap(colvarvalue &x) const;
/// Read the analysis tasks
int parse_analysis(std::string const &conf);
/// Perform analysis tasks
void analyze();
/// Read the value from a collective variable trajectory file
std::istream & read_traj(std::istream &is);
/// Output formatted values to the trajectory file
std::ostream & write_traj(std::ostream &os);
/// Write a label to the trajectory file (comment line)
std::ostream & write_traj_label(std::ostream &os);
/// Read the collective variable from a restart file
std::istream & read_restart(std::istream &is);
/// Write the collective variable to a restart file
std::ostream & write_restart(std::ostream &os);
/// Write output files (if defined, e.g. in analysis mode)
int write_output_files();
protected:
/// Previous value (to calculate velocities during analysis)
colvarvalue x_old;
/// Time series of values and velocities used in correlation
/// functions
std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
/// Time series of values and velocities used in correlation
/// functions (pointers)x
std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;
/// Time series of values and velocities used in running averages
std::list< std::list<colvarvalue> > x_history;
/// Time series of values and velocities used in correlation
/// functions (pointers)x
std::list< std::list<colvarvalue> >::iterator x_history_p;
/// \brief Collective variable with which the correlation is
/// calculated (default: itself)
std::string acf_colvar_name;
/// Length of autocorrelation function (ACF)
size_t acf_length;
/// After how many steps the ACF starts
size_t acf_offset;
/// How many timesteps separate two ACF values
size_t acf_stride;
/// Number of frames for each ACF point
size_t acf_nframes;
/// Normalize the ACF to a maximum value of 1?
bool acf_normalize;
/// ACF values
std::vector<cvm::real> acf;
/// Name of the file to write the ACF
std::string acf_outfile;
/// Type of autocorrelation function (ACF)
enum acf_type_e {
/// Unset type
acf_notset,
/// Velocity ACF, scalar product between v(0) and v(t)
acf_vel,
/// Coordinate ACF, scalar product between x(0) and x(t)
acf_coor,
/// \brief Coordinate ACF, second order Legendre polynomial
/// between x(0) and x(t) (does not work with scalar numbers)
acf_p2coor
};
/// Type of autocorrelation function (ACF)
acf_type_e acf_type;
/// \brief Velocity ACF, scalar product between v(0) and v(t)
int calc_vel_acf(std::list<colvarvalue> &v_history,
colvarvalue const &v);
/// \brief Coordinate ACF, scalar product between x(0) and x(t)
/// (does not work with scalar numbers)
void calc_coor_acf(std::list<colvarvalue> &x_history,
colvarvalue const &x);
/// \brief Coordinate ACF, second order Legendre polynomial between
/// x(0) and x(t) (does not work with scalar numbers)
void calc_p2coor_acf(std::list<colvarvalue> &x_history,
colvarvalue const &x);
/// Calculate the auto-correlation function (ACF)
int calc_acf();
/// Save the ACF to a file
void write_acf(std::ostream &os);
/// Length of running average series
size_t runave_length;
/// Timesteps to skip between two values in the running average series
size_t runave_stride;
/// Name of the file to write the running average
cvm::ofstream runave_os;
/// Current value of the running average
colvarvalue runave;
/// Current value of the square deviation from the running average
cvm::real runave_variance;
/// Calculate the running average and its standard deviation
void calc_runave();
/// If extended Lagrangian active: colvar energies (kinetic and harmonic potential)
cvm::real kinetic_energy;
cvm::real potential_energy;
public:
// collective variable component base class
class cvc;
// currently available collective variable components
// scalar colvar components
class distance;
class distance_z;
class distance_xy;
class distance_inv;
class distance_pairs;
class angle;
class dipole_angle;
class dihedral;
class coordnum;
class selfcoordnum;
class h_bond;
class rmsd;
class orientation_angle;
class orientation_proj;
class tilt;
class spin_angle;
class gyration;
class inertia;
class inertia_z;
class eigenvector;
class alpha_dihedrals;
class alpha_angles;
class dihedPC;
// non-scalar components
class distance_vec;
class distance_dir;
class cartesian;
class orientation;
protected:
/// \brief Array of \link cvc \endlink objects
std::vector<cvc *> cvcs;
/// \brief Flags to enable or disable cvcs at next colvar evaluation
std::vector<bool> cvc_flags;
/// \brief Initialize the sorted list of atom IDs for atoms involved
/// in all cvcs (called when enabling f_cv_collect_gradients)
void build_atom_list(void);
private:
/// Name of scripted function to be used
std::string scripted_function;
/// Current cvc values in the order requested by script
/// when using scriptedFunction
std::vector<const colvarvalue *> sorted_cvc_values;
public:
/// \brief Sorted array of (zero-based) IDs for all atoms involved
std::vector<int> atom_ids;
/// \brief Array of atomic gradients collected from all cvcs
/// with appropriate components, rotations etc.
/// For scalar variables only!
std::vector<cvm::rvector> atomic_gradients;
inline size_t n_components() const {
return cvcs.size();
}
};
inline colvarvalue const & colvar::value() const
{
return x_reported;
}
inline colvarvalue const & colvar::actual_value() const
{
return x;
}
inline colvarvalue const & colvar::velocity() const
{
return v_reported;
}
inline colvarvalue const & colvar::system_force() const
{
return ft_reported;
}
inline void colvar::add_bias_force(colvarvalue const &force)
{
if (cvm::debug()) {
cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
}
fb += force;
}
inline void colvar::reset_bias_force() {
fb.type(value());
fb.reset();
}
#endif
Event Timeline
Log In to Comment