Analyzing source code of Colvars package
In this post we look at some parts of the source code of the Colvars package (module).
Compilation with COLVARS_DEBUG
To get a understanding of the source code, it is helpful to compile the Colvars module with an MD engine with the macro COLVARS_DEBUG set to true and run MD tests. This will instruct Colvars to produce lots of runtime information which is helpful for debugging and analyzing the code. However, for some reason adding the flag -DCOLVARS_DEBUG in compilation does not have an effect in my case. Therefore, I simply modified the first lines of the file colvarmodule.h, defined the macro by force as follows
#ifndef COLVARS_DEBUG
//#define COLVARS_DEBUG false # original line, replaced by the line below
#define COLVARS_DEBUG true # define the macro
and compiled the Colvars module with Gromacs following a previous post. After installation, running an MD simulation of alanine dipeptide with adaptive biasing force (ABF) method produces the following output:
...
colvars: colvarproxy_gromacs, step no. 0
colvars: Updating internal data.
colvars: ----------------------------------------------------------------------
colvars: Collective variables module, step no. 1
colvars: Calculating collective variables.
colvars: Calculating colvar "phi", components 0 through 0.
colvars: Calculating colvar components.
colvars: Colvar component no. 1 within colvar "phi" has value -1.35499242515650e+02.
colvars: Calculating gradients of colvar "phi".
colvars: Done calculating gradients of colvar "phi".
colvars: Done calculating colvar "phi".
colvars: Calculating colvar "phi"'s properties.
colvars: Colvar "phi" has value -1.35499242515650e+02.
colvars: Done calculating colvar "phi"'s properties.
colvars: Calculating colvar "psi", components 0 through 0.
colvars: Calculating colvar components.
colvars: Colvar component no. 1 within colvar "psi" has value 1.65590822749484e+02.
colvars: Calculating gradients of colvar "psi".
colvars: Done calculating gradients of colvar "psi".
colvars: Done calculating colvar "psi".
colvars: Calculating colvar "psi"'s properties.
colvars: Colvar "psi" has value 1.65590822749484e+02.
colvars: Done calculating colvar "psi"'s properties.
colvars: Updating collective variable biases.
colvars: Updating ABF bias abf1
colvars: ABF System force calc: cv 0 fs 0 = ft 0 - fa 0
colvars: ABF System force calc: cv 1 fs 0 = ft 0 - fa 0
colvars: ABF System force calc: cv 0 fs 0 = ft 0 - fa 0
colvars: ABF System force calc: cv 1 fs 0 = ft 0 - fa 0
colvars: Collecting forces from all biases.
colvars: Communicating a force to colvar "phi".
colvars: Adding biasing force 0 to colvar "phi".
colvars: Communicating a force to colvar "psi".
colvars: Adding biasing force 0 to colvar "psi".
colvars: Adding total bias energy: 0
colvars: Updating the internal degrees of freedom of colvars (if they have any).
colvars: Updating colvar "phi".
colvars: Updating extended-Lagrangian degree of freedom.
colvars: Done updating colvar "phi".
colvars: Updating colvar "psi".
colvars: Updating extended-Lagrangian degree of freedom.
colvars: Done updating colvar "psi".
colvars: Adding total colvar energy: 0.0302533
colvars: Communicating forces from the colvars to the atoms.
colvars: Communicating forces from colvar "phi".
colvars: Force to be applied: 0.0341128
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Done communicating forces from colvar "phi".
colvars: Communicating forces from colvar "psi".
colvars: Force to be applied: -0.0638968
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Communicating a colvar force from atom group to the MD engine.
colvars: Done communicating forces from colvar "psi".
colvars: colvarmodule::analyze(), step = 1.
colvars: colvarmodule::write_traj_files()
colvars: Using colvarproxy_io::output_stream()
colvars: colvarmodule::end_of_step(), step = 1.
colvars: End of step for colvar "phi".
colvars: End of step for colvar "psi".
colvars: ----------------------------------------------------------------------
...
In the following, we analyze the source code of Colvars package and identify functions which generate the above information.
Level 1. classes colvarproxy_*
These are classes defined to interface between Colvars module and MD engine.
For example, one can find in the header file colvarproxy.h the following code:
/// \brief Container of atomic data for processing by Colvars
class colvarproxy_atoms {
...
/// Request that this force is applied to the given atom
/// \param index Internal index in the Colvars arrays
/// \param new_force Force to add
inline void apply_atom_force(int index, cvm::rvector const &new_force)
{
atoms_new_colvar_forces[index] += new_force;
}
...
protected:
/// \brief Array of 0-based integers used to uniquely associate atoms
/// within the host program
std::vector<int> atoms_ids;
/// \brief Keep track of how many times each atom is used by a separate colvar object
std::vector<size_t> atoms_refcount;
/// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_masses;
/// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
std::vector<cvm::real> atoms_charges;
/// \brief Current three-dimensional positions of the atoms
std::vector<cvm::rvector> atoms_positions;
/// \brief Most recent total forces on each atom
std::vector<cvm::rvector> atoms_total_forces;
/// \brief Forces applied from colvars, to be communicated to the MD integrator
std::vector<cvm::rvector> atoms_new_colvar_forces;
...
};
...
/// \brief Container of atom group data (allow collection of aggregated atomic
/// data)
class colvarproxy_atom_groups {
...
/// Request that this force is applied to the given atom group
inline void apply_atom_group_force(int index, cvm::rvector const &new_force)
{
atom_groups_new_colvar_forces[index] += new_force;
}
...
protected:
/// \brief Array of 0-based integers used to uniquely associate atom groups
/// within the host program
std::vector<int> atom_groups_ids;
/// \brief Keep track of how many times each group is used by a separate cvc
std::vector<size_t> atom_groups_refcount;
/// \brief Total masses of the atom groups
std::vector<cvm::real> atom_groups_masses;
...
};
/// Interface between Colvars and MD engine (GROMACS, LAMMPS, NAMD, VMD...)
///
/// This is the base class: each engine is supported by a derived class.
class colvarproxy
: public colvarproxy_system,
public colvarproxy_atoms,
public colvarproxy_atom_groups,
public colvarproxy_volmaps,
public colvarproxy_smp,
public colvarproxy_replicas,
public colvarproxy_script,
public colvarproxy_tcl,
public colvarproxy_io
{
...
};
In an MD simulation, Colvars module communicates with MD engine through derived classes of class colvarproxy. When the MD engine is Gromacs, the derived interfacing class is colvarproxy_gromacs, which is defined in the hearder file colvarproxy_gromacs.h:
/// \brief Communication between colvars and Gromacs (implementation of
/// \link colvarproxy \endlink)
class colvarproxy_gromacs : public colvarproxy, public gmx::IForceProvider {
...
virtual void calculateForces(const gmx::ForceProviderInput &forceProviderInput,
gmx::ForceProviderOutput *forceProviderOutput);
...
};
The implementation is in colvarproxy_gromacs.cpp. In particular, the Colvars module is invoked in the following function.
void colvarproxy_gromacs::calculateForces(
const gmx::ForceProviderInput &forceProviderInput,
gmx::ForceProviderOutput *forceProviderOutput)
{
const t_commrec *cr = &(forceProviderInput.cr_);
// Local atom coords
const gmx::ArrayRef<const gmx::RVec> x = forceProviderInput.x_;
// Local atom coords (coerced into into old gmx type)
const rvec *x_pointer = &(x.data()->as_vec());
// Eventually there needs to be an interface to update local data upon neighbor search
// We could check if by chance all atoms are in one node, and skip communication
communicate_group_positions(cr, x_colvars_unwrapped, xa_shifts, xa_eshifts,
gmx_bNS, x_pointer, n_colvars_atoms, nat_loc,
ind_loc, xa_ind, xa_old_whole, gmx_box);
// Communicate_group_positions takes care of removing shifts (unwrapping)
// in single node jobs, communicate_group_positions() is efficient and adds no overhead
if (MASTER(cr))
{
// On non-master nodes, jump directly to applying the forces
// Zero the forces on the atoms, so that they can be accumulated by the colvars.
for (size_t i = 0; i < atoms_new_colvar_forces.size(); i++) {
atoms_new_colvar_forces[i].x = atoms_new_colvar_forces[i].y = atoms_new_colvar_forces[i].z = 0.0;
}
// Get the atom positions from the Gromacs array.
for (size_t i = 0; i < atoms_ids.size(); i++) {
atoms_positions[i] = cvm::rvector(x_colvars_unwrapped[i][0], x_colvars_unwrapped[i][1], x_colvars_unwrapped[i][2]);
}
bias_energy = 0.0;
// Call the collective variable module to fill atoms_new_colvar_forces
if (colvars->calc() != COLVARS_OK) {
cvm::error("Error calling colvars->calc()\n");
}
// Copy the forces to C array for broadcasting
for (int i = 0; i < n_colvars_atoms; i++)
{
f_colvars[i][0] = atoms_new_colvar_forces[i].x;
f_colvars[i][1] = atoms_new_colvar_forces[i].y;
f_colvars[i][2] = atoms_new_colvar_forces[i].z;
}
forceProviderOutput->enerd_.term[F_COM_PULL] += bias_energy;
} // master node
//Broadcast the forces to all the nodes
if (PAR(cr))
{
nblock_bc(cr->mpi_comm_mygroup, n_colvars_atoms, f_colvars);
}
const gmx::ArrayRef<gmx::RVec> &f_out = forceProviderOutput->forceWithVirial_.force_;
matrix local_colvars_virial = { { 0 } };
// Pass the applied forces back to GROMACS
for (int i = 0; i < n_colvars_atoms; i++)
{
int i_global = ind[i];
// check if this is a local atom and find out locndx
if (PAR(cr)) {
const int *locndx = cr->dd->ga2la->findHome(i_global);
if (locndx) {
f_out[*locndx] += f_colvars[i];
add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]);
}
// Do nothing if atom is not local
} else { // Non MPI-parallel
f_out[i_global] += f_colvars[i];
add_virial_term(local_colvars_virial, f_colvars[i], x_colvars_unwrapped[i]);
}
}
forceProviderOutput->forceWithVirial_.addVirialContribution(local_colvars_virial);
return;
}
From the function above, one can see that the Colvars module first retrieves position data of the system from Gromacs before it starts calculation, by calling the following function defined in src/gromacs/mdlib/groupcoord.cpp (under the home directory of Gromacs package).
/* Assemble the positions of the group such that every node has all of them.
* The atom indices are retrieved from anrs_loc[0..nr_loc]
* Note that coll_ind[i] = i is needed in the serial case */
extern void communicate_group_positions(const t_commrec* cr, /* Pointer to MPI communication data */
rvec* xcoll, /* Collective array of positions */
ivec* shifts, /* Collective array of shifts for xcoll (can be NULL) */
ivec* extra_shifts, /* (optional) Extra shifts since last time step */
const gmx_bool bNS, /* (optional) NS step, the shifts have changed */
const rvec* x_loc, /* Local positions on this node */
const int nr, /* Total number of atoms in the group */
const int nr_loc, /* Local number of atoms in the group */
const int* anrs_loc, /* Local atom numbers */
const int* coll_ind, /* Collective index */
rvec* xcoll_old, /* (optional) Positions from the last time
step, used to make group whole */
const matrix box) /* (optional) The box */
{
int i;
/* Zero out the groups' global position array */
clear_rvecs(nr, xcoll);
/* Put the local positions that this node has into the right place of
* the collective array. Note that in the serial case, coll_ind[i] = i */
for (i = 0; i < nr_loc; i++)
{
copy_rvec(x_loc[anrs_loc[i]], xcoll[coll_ind[i]]);
}
if (PAR(cr))
{
/* Add the arrays from all nodes together */
gmx_sum(nr * 3, xcoll[0], cr);
}
/* Now we have all the positions of the group in the xcoll array present on all
* nodes.
*
* The rest of the code is for making the group whole again in case atoms changed
* their PBC representation / crossed a box boundary. We only do that if the
* shifts array is allocated. */
if (nullptr != shifts)
{
/* To make the group whole, start with a whole group and each
* step move the assembled positions at closest distance to the positions
* from the last step. First shift the positions with the saved shift
* vectors (these are 0 when this routine is called for the first time!) */
shift_positions_group(box, xcoll, shifts, nr);
/* Now check if some shifts changed since the last step.
* This only needs to be done when the shifts are expected to have changed,
* i.e. after neighbor searching */
if (bNS)
{
get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
/* Shift with the additional shifts such that we get a whole group now */
shift_positions_group(box, xcoll, extra_shifts, nr);
/* Add the shift vectors together for the next time step */
for (i = 0; i < nr; i++)
{
shifts[i][XX] += extra_shifts[i][XX];
shifts[i][YY] += extra_shifts[i][YY];
shifts[i][ZZ] += extra_shifts[i][ZZ];
}
/* Store current correctly-shifted positions for comparison in the next NS time step */
for (i = 0; i < nr; i++)
{
copy_rvec(xcoll[i], xcoll_old[i]);
}
}
}
}
Level 2. class colvarmodule
This class is the main class of Colvars module. It invokes the executions of various computational tasks related to collective varialbes (see the file colvarmodule.cpp). In particular, the function below is called by the class colvarproxy_gromacs above.
2.1. Main function to invoke various computations
int colvarmodule::calc()
{
int error_code = COLVARS_OK;
if (cvm::debug()) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables module, step no. "+
cvm::to_str(cvm::step_absolute())+"\n");
}
error_code |= calc_colvars();
error_code |= calc_biases();
error_code |= update_colvar_forces();
error_code |= analyze();
...
}
2.2. Function to compute values and gradients of CVs
int colvarmodule::calc_colvars()
{
if (cvm::debug())
cvm::log("Calculating collective variables.\n");
...
for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
error_code |= (*cvi)->calc();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
...
}
This function calls the function int colvar::calc() to perform actual computations. See Section 3.2.1.1 below.
2.3. Function to updates the FE gradient, compute and apply the biasing forces
int colvarmodule::calc_biases()
{
...
if (cvm::debug() && num_biases())
cvm::log("Updating collective variable biases.\n");
...
for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
error_code |= (*bi)->update();
if (cvm::get_error()) {
return error_code;
}
}
...
}
This function calls int colvarbias_abf::update() (or functions of other biasing algorithms) to perform the actual computations. See Section 3.1.1 below.
2.4. Function to update CV forces from biasing methods and apply forces to atoms
int colvarmodule::update_colvar_forces()
{
if (cvm::debug() && num_biases())
cvm::log("Collecting forces from all biases.\n");
...
for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
error_code |= (*bi)->communicate_forces();
}
...
for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) {
// Inactive colvars will only reset their forces and return 0 energy
total_colvar_energy += (*cvi)->update_forces_energy();
}
...
for (cvi = variables_active()->begin(); cvi != variables_active()->end(); cvi++) {
if ((*cvi)->is_enabled(colvardeps::f_cv_gradient)) {
(*cvi)->communicate_forces();
if (cvm::get_error()) {
return COLVARS_ERROR;
}
}
}
}
This function calls the function int colvarbias::communicate_forces() to update the CV forces (Section 3.1.4) and the function void colvar::communicate_forces() to apply the forces to atoms (Section 3.2.4), respectively.
2.5. atom and atom groups
Besides various computational tasks, colvarmodule class also declares classes for atom and atom group in the header file colvarmodule.h:
class colvarmodule {
...
// allow these classes to access protected data
class atom;
class atom_group;
friend class atom;
friend class atom_group;
...
}
These two classes are defined in colvaratoms.h and colvaratoms.cpp.
In particular, colvar forces are communicated from atom group to MD engine by the following two functions (See the function colvar::dihedral::apply_force in the Level 4 below).
void cvm::atom_group::apply_colvar_force(cvm::real const &force)
{
if (cvm::debug()) {
log("Communicating a colvar force from atom group to the MD engine.\n");
}
...
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
return;
}
if (is_enabled(f_ag_rotate)) {
// rotate forces back to the original frame
cvm::rotation const rot_inv = rot.inverse();
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->apply_force(rot_inv.rotate(force * ai->grad));
}
} else {
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->apply_force(force * ai->grad);
}
}
...
}
void cvm::atom_group::apply_force(cvm::rvector const &force)
{
if (cvm::debug()) {
log("Communicating a colvar force from atom group to the MD engine.\n");
}
...
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force);
return;
}
if (is_enabled(f_ag_rotate)) {
cvm::rotation const rot_inv = rot.inverse();
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->apply_force(rot_inv.rotate((ai->mass/total_mass) * force));
}
} else {
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
ai->apply_force((ai->mass/total_mass) * force);
}
}
}
These functions call the following one (defined in the header file colvaratoms.h) to communicate forces to each atom:
class colvarmodule::atom {
...
/// \brief Apply a force to the atom
///
/// Note: the force is not applied instantly, but will be used later
/// by the MD integrator (the colvars module does not integrate
/// equations of motion.
///
/// Multiple calls to this function by either the same
/// \link atom \endlink object or different objects with identical
/// \link id \endlink will all be added together.
inline void apply_force(cvm::rvector const &new_force) const
{
(cvm::proxy)->apply_atom_force(index, new_force);
}
};
This function eventually calls the function cvm::proxy::apply_atom_force defined in class colvarproxy_atom to communciate with MD engine (e.g. Gromacs). See the beginning of Level 1 above.
Level 3. class colvarbias_abf and class colvar
Level 3.1. class colvarbias_abf
colvarbias is the base class for biasing algorihms. The ABF algorithm is implemented in the derived class colvarbias_abf.
3.1.1. Function to update the FE gradient and compute CV forces
This function is invoked by the function int colvarmodule::calc_biases() in Section 2.3 above.
int colvarbias_abf::update()
{
if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
...
if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
if (update_bias) {
if (samples->index_ok(force_bin)) {
// Only if requested and within bounds of the grid...
for (i = 0; i < num_variables(); i++) {
// get total forces (lagging by 1 timestep) from colvars
// and subtract previous ABF force if necessary
update_system_force(i);
}
gradients->acc_force(force_bin, system_force);
}
}
...
}
...
// Compute and apply the new bias, if applicable
if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
...
// Factor that ensures smooth introduction of the force
if ( count < full_samples ) {
fact = (count < min_samples) ? 0.0 :
(cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
}
...
if ( pabf_freq ) {
// In projected ABF, the force is the PMF gradient estimate
pmf->vector_gradient_finite_diff(bin, grad);
} else {
// Normal ABF
gradients->vector_value(bin, grad);
}
if ( fact != 0.0 ) {
if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) {
// Enforce a zero-mean bias on periodic, 1D coordinates
// in other words: boundary condition is that the biasing potential is periodic
// This is enforced naturally if using integrated PMF
colvar_forces[0].real_value = fact * (grad[0] - gradients->average ());
} else {
for (i = 0; i < num_variables(); i++) {
// subtracting the mean force (opposite of the FE gradient) means adding the gradient
colvar_forces[i].real_value = fact * grad[i];
}
}
}
...
}
}
This function calls the functions in Section 3.1.2. and Section 3.1.3 below.
3.1.2. Function to update system force from colvar
(defined in colvarbias_abf.h)
inline int update_system_force(size_t i)
{
if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) {
// this colvar is already subtracting the ABF force
system_force[i] = colvars[i]->total_force().real_value;
} else {
system_force[i] = colvars[i]->total_force().real_value
- colvar_forces[i].real_value;
// If hideJacobian is active then total_force has an extra term of -fj
// which is the Jacobian-compensating force at the colvar level
}
if (cvm::debug())
cvm::log("ABF System force calc: cv " + cvm::to_str(i) +
" fs " + cvm::to_str(system_force[i]) +
" = ft " + cvm::to_str(colvars[i]->total_force().real_value) +
" - fa " + cvm::to_str(colvar_forces[i].real_value));
return COLVARS_OK;
}
3.1.3. Function to accumulate force on grid mesh
(defined in colvargrid.h)
/// \brief Accumulate the gradient based on the force (i.e. sums the
/// opposite of the force)
inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
for (size_t imult = 0; imult < mult; imult++) {
data[address(ix) + imult] -= forces[imult];
}
if (samples)
samples->incr_count(ix);
}
3.1.4. Function to communicate CV forces (computed by biasing methods) to class colvar
int colvarbias::communicate_forces()
{
...
for (i = 0; i < num_variables(); i++) {
if (cvm::debug()) {
cvm::log("Communicating a force to colvar \""+
variables(i)->name+"\".\n");
}
// Impulse-style multiple timestep
// Note that biases with different values of time_step_factor
// may send forces to the same colvar
// which is why rescaling has to happen now: the colvar is not
// aware of this bias' time_step_factor
if (is_enabled(f_cvb_bypass_ext_lagrangian)) {
variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i] * biasing_force_factor);
} else {
variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i] * biasing_force_factor);
}
previous_colvar_forces[i] = colvar_forces[i];
}
}
This function is invoked by the function int colvarmodule::update_colvar_forces() in Section 2.4 above. It in turn calls the function inline void colvar::add_bias_force(colvarvalue const &force) to add bias forces to CVs (in Section 3.2.3).
Level 3.2. class colvar
This class implements various calculations related to CVs.
3.2.1. Functions to perform CV-related computations
3.2.1.1. Main function
The main function is invoked by the function int colvarmodule::calc_colvars() in Section 2.2.
// Default schedule (everything is serialized)
int colvar::calc()
{
// Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
int error_code = COLVARS_OK;
if (is_enabled(f_cv_active)) {
error_code |= update_cvc_flags();
if (error_code != COLVARS_OK) return error_code;
error_code |= calc_cvcs();
if (error_code != COLVARS_OK) return error_code;
error_code |= collect_cvc_data();
}
return error_code;
}
3.2.1.2. CV computations
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\", components "+
cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
...
if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
// Use Jacobian derivative from previous timestep
error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
}
// atom coordinates are updated by the next line
error_code |= calc_cvc_values(first_cvc, num_cvcs);
error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
if (proxy->total_forces_same_step()){
// Use Jacobian derivative from this timestep
error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
}
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\".\n");
return error_code;
}
int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs)
{
size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
size_t i, cvc_count;
// calculate the value of the colvar
if (cvm::debug())
cvm::log("Calculating colvar components.\n");
// First, calculate component values
cvm::increase_depth();
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->is_enabled()) continue;
cvc_count++;
(cvcs[i])->read_data();
(cvcs[i])->calc_value();
if (cvm::debug())
cvm::log("Colvar component no. "+cvm::to_str(i+1)+
" within colvar \""+this->name+"\" has value "+
cvm::to_str((cvcs[i])->value(),
cvm::cv_width, cvm::cv_prec)+".\n");
}
cvm::decrease_depth();
return COLVARS_OK;
}
int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
{
size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
size_t i, cvc_count;
if (cvm::debug())
cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
// calculate the gradients of each component
cvm::increase_depth();
for (i = first_cvc, cvc_count = 0;
(i < cvcs.size()) && (cvc_count < cvc_max_count);
i++) {
if (!cvcs[i]->is_enabled()) continue;
cvc_count++;
if ((cvcs[i])->is_enabled(f_cvc_gradient)) {
(cvcs[i])->calc_gradients();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
(cvcs[i])->calc_fit_gradients();
if ((cvcs[i])->is_enabled(f_cvc_debug_gradient))
(cvcs[i])->debug_gradients();
}
cvm::decrease_depth();
if (cvm::debug())
cvm::log("Done calculating gradients of colvar \""+this->name+"\".\n");
}
return COLVARS_OK;
}
These functions call functions of Colvar components (see Level 4 below) to perform actual computations.
3.2.1.3. Collect data
int colvar::collect_cvc_data()
{
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
// Total force depends on Jacobian derivative from previous timestep
// collect_cvc_total_forces() uses the previous value of jd
error_code |= collect_cvc_total_forces();
}
error_code |= collect_cvc_values();
error_code |= collect_cvc_gradients();
error_code |= collect_cvc_Jacobians();
if (proxy->total_forces_same_step()){
// Use Jacobian derivative from this timestep
error_code |= collect_cvc_total_forces();
}
error_code |= calc_colvar_properties();
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
return error_code;
}
3.2.2. Update force and energy
cvm::real colvar::update_forces_energy()
{
if (cvm::debug())
cvm::log("Updating colvar \""+this->name+"\".\n");
// add the biases' force, which at this point should already have
// been summed over each bias using this colvar
// fb is already multiplied by the relevant time step factor for each bias
f += fb;
...
// At this point f is the force f from external biases that will be applied to the
// extended variable if there is one
if (is_enabled(f_cv_extended_Lagrangian) && cvm::proxy->simulation_running()) {
update_extended_Lagrangian();
}
...
}
This function is called by int colvarmodule::update_colvar_forces() in Section 2.4. When the flag extendedLagrangian is true for a colvar, this function in turns call void colvar::update_extended_Lagrangian() to update the extended degree of freedom. The corresponding code in colvar.cpp is worth reading.
3.2.3. Add bias forces
class colvar : public colvarparse, public colvardeps {
...
/// Add to the total force from biases
void add_bias_force(colvarvalue const &force);
/// Apply a force to the actual value (only meaningful with extended Lagrangian)
void add_bias_force_actual_value(colvarvalue const &force);
...
}
inline void colvar::add_bias_force(colvarvalue const &force)
{
check_enabled(f_cv_gradient,
std::string("applying a force to the variable \""+name+"\""));
if (cvm::debug()) {
cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
}
fb += force;
}
inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
{
if (cvm::debug()) {
cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
}
fb_actual += force;
}
These functions are called by int colvarbias::communicate_forces() in Section 3.1.4.
3.2.4. Communicate forces from a CV to atoms
void colvar::communicate_forces()
{
if (cvm::debug()) {
cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
cvm::log("Force to be applied: " + cvm::to_str(f) + "\n");
}
...
for (i = 0; i < cvcs.size(); i++) {
if (!cvcs[i]->is_enabled()) continue;
(cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
cvm::real((cvcs[i])->sup_np) *
(cvm::integer_power((cvcs[i])->value().real_value,
(cvcs[i])->sup_np-1)) );
}
This function is called by int colvarmodule::update_colvar_forces() in Section 2.4. It invokes the function of Colvars component (in Level 4 below) to apply forces to atoms.
Level 4. class colvar::cvc
This class describes a Colvar component. An object of this class (or an object of a cvc-derived class) implements the calculation of a collective variable and its gradients. Many concrete colvar components, such as distances, angles, and dihedral angles, are derived from this class. They are defined in the header file colvarcomp.h and in the cpp files.
For example, below is the code in the header file for an angle component.
/// \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;
...
}
The function called by the colvar::communicate_forces() (see Section 3.2.4. above) is implemented in colvarcomp_angles.cpp as follows:
void colvar::angle::apply_force(colvarvalue const &force)
{
if (!group1->noforce)
group1->apply_colvar_force(force.real_value);
if (!group2->noforce)
group2->apply_colvar_force(force.real_value);
if (!group3->noforce)
group3->apply_colvar_force(force.real_value);
}
The two functions below are called by the functions int colvar::calc_cvc_values(int first_cvc, size_t num_cvcs) and int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs) in Section 3.2.1.2.
void colvar::angle::calc_value()
{
cvm::atom_pos const g1_pos = group1->center_of_mass();
cvm::atom_pos const g2_pos = group2->center_of_mass();
cvm::atom_pos const g3_pos = group3->center_of_mass();
r21 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g1_pos) :
g1_pos - g2_pos;
r21l = r21.norm();
r23 = is_enabled(f_cvc_pbc_minimum_image) ?
cvm::position_distance(g2_pos, g3_pos) :
g3_pos - g2_pos;
r23l = r23.norm();
cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
x.real_value = (180.0/PI) * cvm::acos(cos_theta);
}
void colvar::angle::calc_gradients()
{
cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
cvm::real const dxdcos = -1.0 / cvm::sqrt(1.0 - cos_theta*cos_theta);
dxdr1 = (180.0/PI) * dxdcos *
(1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l );
dxdr3 = (180.0/PI) * dxdcos *
(1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
group1->set_weighted_gradient(dxdr1);
group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
group3->set_weighted_gradient(dxdr3);
}