Implementation of OPES in Colvars

[HanatoK]

This PR intends to implement the OPES method in Colvars as requested by Chris Chipot. Basically I just followed all the logic in OPESmetad.cpp of the PLUMED implementation, with some Colvars specific modifications.

Things still need to be done (Done):

• Add a basic test
• Add a test for restart
• Add documentation
• Colvars does not have the feature of "chaining two biases" or "reusing some output components of a bias in another bias". The PMF can be otained on-the-fly by chaining OPES_METAD and REWEIGHTING_BIAS in PLUMED, but that is not possible in Colvars, so I will have to implement a weighted histogram inside the bias. This implementation, however, still has two limitations: (i) one cannot apply bias on $(\xi_1, \xi_2)$ and then estimate the PMF along $\xi_3$ by reweighting, and (ii) one cannot estimate the PMFs along $\xi_1$ and $\xi_2$ simultaneously while applying bias on $(\xi_1, \xi_2)$.

Advanced features (and any of their combinations):

• Neighbor list
• Multiple walker OPES
• Adaptive sigma
• Exploration mode

Limitations:

• Due to the same limitation above, the OPES implementation in Colvars cannot take the biasing energies of other biases into account;
• PLUMED does not calculate actions in parallel and therefore allows the OpenMP threads to be used to sum the kernels in parallel. In contrary, Colvars computes the biases themselves on-the-fly, but it is not clear to me if the threads can be used in a single bias to parallelize some of the computations. Although I have written the code, it still needs more testing.

[HanatoK]

Hi @giacomofiorin and @jhenin! It seems that in the latex document the default value cannot be an equation. Otherwise the HTML files are not generated. For example, I want to add the following key description:

\item %
  \labelkey{opes_metad|kernel_cutoff}
  \keydef
    {kernel{\textunderscore}cutoff}{%
    \texttt{opes{\textunderscore}metad}}{%
    Truncate kernels at this distance, in units of $\sigma$}{%
    positive decimal}{%
    $\sqrt{2*\mathrm{barrier}/k_\mathrm{B}T/(1-1/\gamma)}$}{%
    If the value of $(\boldsymbol{\xi}-\boldsymbol{\xi}')^T \Sigma^{-1}(\boldsymbol{\xi}-\boldsymbol{\xi}')$ is larger than the square of this value, then $G(\boldsymbol{\xi}, \boldsymbol{\xi}')$ is evaluated as zero.}

What is the proper way to do it to generate the HTML document?

Update:
Please ignore the question. It seems that the issue is caused by my mistake of mixing \key and \keydef.

[invemichele]

Thanks for bringing OPES to Colvars!

[HanatoK]

Thanks for bringing OPES to Colvars!

Nice to see you! Your idea of OPES is really amazing!

[giacomofiorin]

Update: Please ignore the question. It seems that the issue is caused by my mistake of mixing \key and \keydef.

We should probably rename them to something clearer.

[giacomofiorin]

This is a quick round looking mostly at the documentation. The added source file is very large, so I'll take a look at the differences between PLUMED's version and this one.

One thing that really stands out is that the user is presented with a lot of restrictions, which for the most part come from the fact that biases are distributed among SMP threads. However, there isn't a huge performance gain from distributing biases. I think it would be okay to run all biases in sequence whenever a special dependency is present, as we do already in several cases:

colvars/src/colvarmodule.cpp

Lines 1016 to 1042 in a587282

	if (proxy->check_smp_enabled() == COLVARS_OK && !biases_need_main_thread) {
	if (use_scripted_forces && !scripting_after_biases) {
	// calculate biases and scripted forces in parallel
	error_code |= proxy->smp_biases_script_loop();
	} else {
	// calculate biases in parallel
	error_code |= proxy->smp_biases_loop();
	}
	} else {
	if (use_scripted_forces && !scripting_after_biases) {
	error_code |= calc_scripted_forces();
	}
	// Straight loop over biases on a single thread
	cvm::increase_depth();
	for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
	error_code |= (*bi)->update();
	if (cvm::get_error()) {
	cvm::decrease_depth();
	return error_code;
	}
	}
	cvm::decrease_depth();
	}

I'll comment on #712: would implementing that change allow to remove the restrictions that this OPES implementation has compared to the PLUMED one, which people are already familiar with?

[giacomofiorin]

Is multiple-walkers supported without MPI/Charm++?

[HanatoK]

No. Multiple-walkers OPES requires MPI or Charm++.

[giacomofiorin]

If FILE is not needed because Colvars restarts from a single checkpoint file, then this could be maybe clarified somewhere?

[HanatoK]

Done. (see 3dec310)

[giacomofiorin]

This keyword is listed as mandatory/compulsory in PLUMED, so maybe it shouldn't have a default value?

[HanatoK]

No. Compulsory keywords in PLUMED could have default values (see https://www.plumed.org/doc-master/user-doc/html/_o_p_e_s__m_e_t_a_d__e_x_p_l_o_r_e.html):

COMPRESSION_THRESHOLD 	( default=1 ) merge kernels if closer than this threshold, in units of sigma 

[giacomofiorin]

Thanks!

[HanatoK]

This is a quick round looking mostly at the documentation. The added source file is very large, so I'll take a look at the differences between PLUMED's version and this one.

One thing that really stands out is that the user is presented with a lot of restrictions, which for the most part come from the fact that biases are distributed among SMP threads. However, there isn't a huge performance gain from distributing biases. I think it would be okay to run all biases in sequence whenever a special dependency is present, as we do already in several cases:

colvars/src/colvarmodule.cpp

Lines 1016 to 1042 in a587282

	if (proxy->check_smp_enabled() == COLVARS_OK && !biases_need_main_thread) {
	if (use_scripted_forces && !scripting_after_biases) {
	// calculate biases and scripted forces in parallel
	error_code |= proxy->smp_biases_script_loop();
	} else {
	// calculate biases in parallel
	error_code |= proxy->smp_biases_loop();
	}
	} else {
	if (use_scripted_forces && !scripting_after_biases) {
	error_code |= calc_scripted_forces();
	}
	// Straight loop over biases on a single thread
	cvm::increase_depth();
	for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) {
	error_code |= (*bi)->update();
	if (cvm::get_error()) {
	cvm::decrease_depth();
	return error_code;
	}
	}
	cvm::decrease_depth();
	}

I'll comment on #712: would implementing that change allow to remove the restrictions that this OPES implementation has compared to the PLUMED one, which people are already familiar with?

Whether the code runs in parallel by OpenMP/ckloop or in serial only affects the performance, and switching between serial and parallel implementations is designed to be almost transparent, although I haven't tested the OpenMP/ckloop code paths yet. I am not sure how you plan to further parallelize Colvars, although I personally incline to a more fine-grained approach of parallelization (using SMP threads for calculating rotations, summing hills and possibly building neighbor list). That being said, this PR only supports the serial code path for the time being, but once you make a decision to switch to a more fine-grained SMP parallelization, I could quickly test the parallelized code path.

[giacomofiorin]

Whether the code runs in parallel by OpenMP/ckloop or in serial only affects the performance, and switching between serial and parallel implementations is designed to be almost transparent,

Performance aside, I was referring specifically to the ability to "chain two biases", which should be straightforward to achieve as long as the order of computation is respected. Would that ability allow a more complete (and thus easier to use) OPES implementation?

[HanatoK]

Whether the code runs in parallel by OpenMP/ckloop or in serial only affects the performance, and switching between serial and parallel implementations is designed to be almost transparent,

Performance aside, I was referring specifically to the ability to "chain two biases", which should be straightforward to achieve as long as the order of computation is respected. Would that ability allow a more complete (and thus easier to use) OPES implementation?

I think the current OPES implementation should be OK for common scenarios. If one wants to apply the biasing potential on $\xi_1$ while expects a PMF along $\xi_2$ or just a free-energy difference by the on-the-fly reweighting, then it is not possible in Colvars, but in PLUMED that is possible by chaining OPES_METAD and REWEIGHT_BIAS. However, I consider this use case is less common, so I have covered the most common use case (biasing $\xi_1$ and expecting a PMF along $\xi_1$) by incorporating an on-the-fly reweighting histogram into the PR. It is also possible to just increase the printing frequency of the trajectory and perform the reweighting after the simulation, if the disk space permits.

As for chaining two biases, I think it would be really useful, not only for OPES, but also for metadynamics. The metadynamics in Colvars currently does not support the calculation of $c(t)$ as described in https://pubs.acs.org/doi/10.1021/jp504920s. The $c(t)$ can be used for reweighting in metadynamics just like the bias energy in OPES. By reweighting one can get rid of reconstructing the full multidimensional PMF and just estimate the unbiased probabilities between two states. Also the free energy surface from reweighting could be converged faster than summing the hills (mentioned in the OPES paper https://dx.doi.org/10.1021/acs.jpclett.0c00497). In PLUMED it is possible to chain "REWEIGHT_METAD" after "METAD" but in Colvars it is currently not possible.