Re: Colvars: implement a CV which requires computing minimal RMSDs with respect to multiple different reference frames

From: Giacomo Fiorin (giacomo.fiorin_at_gmail.com)
Date: Mon Mar 11 2019 - 10:39:31 CDT

Hi Haochuan,

On Mon, Mar 11, 2019 at 9:51 AM <yjcoshc_at_gmail.com> wrote:

> Dear Jerome,
>
> It seems difficult to do something like accepting a vector of RMSDs and
> then applying force on them since RMSD is already a colvar class, and I
> know developers are trying to avoid something like "colvar of colvar".
>
No developers are trying to avoid this as far as I'm aware of.

At the C++ level, "colvar" is the class whose job is controlling how a
variable is *used* (i.e. handling input/output, receive forces from biases,
etc). Instead, colvar::cvc and its derived classes (AKA the components)
are responsible for computing functions of the atomic coordinates. To
compute a function of this function, you are already well aware of the
available methods (deriving a new C++ class from the original one, using a
Tcl scripted function, or a Lepton custom function). Especially if you
derive a new C++ class (e.g. from colvar::rmsd to colvar::multi_rmsd?)
there are hardly any restrictions to what function you can implement.

In your specific example, Jérôme pointed out the easiest solution: a Tcl
proc that computes the minimum of the given arguments (the values of each
RMSD) using a smoothing function of your choice (I presume that you're not
applying forces discontinuously). You'll have, of course, more freedom
with a new C++ class, which you can also document separately from the
original class, and include any references to its use/development in a
paper.

> Today I have managed to make the previous code in my email works by
> modifying colvarparse.cpp like this:
>
Modifications to the parser are tricky, but can be generally useful
provided that your changes are consistent with the rest of the
functionality beyond your use case. Specifically, please ensure that you
are making changes onto the current master (see header at the top of each
Colvars source file).

Giacomo

> diff --git a/src/colvarparse.cpp b/src/colvarparse.cpp
> index 999677e..4d3e8fe 100644
> --- a/src/colvarparse.cpp
> +++ b/src/colvarparse.cpp
> @@ -3,6 +3,7 @@
>
> #include <sstream>
> #include <iostream>
> +#include <algorithm>
>
> #include "colvarmodule.h"
> #include "colvarvalue.h"
> @@ -481,6 +482,10 @@ void colvarparse::strip_values(std::string &conf)
> size_t offset = 0;
> data_begin_pos.sort();
> data_end_pos.sort();
> + std::list<size_t>::iterator data_begin_pos_last =
> std::unique(data_begin_pos.begin(), data_begin_pos.end());
> + data_begin_pos.erase(data_begin_pos_last, data_begin_pos.end());
> + std::list<size_t>::iterator data_end_pos_last =
> std::unique(data_end_pos.begin(), data_end_pos.end());
> + data_end_pos.erase(data_end_pos_last, data_end_pos.end());
>
> std::list<size_t>::iterator data_begin = data_begin_pos.begin();
> std::list<size_t>::iterator data_end = data_end_pos.begin();
>
> The patch just deduplicates the indices when stripping conf strings and
> allows to parse the same keyword repeatedly.
>
> Thanks,
>
> Haochuan Chen
>
>
> 在 2019/3/11 下午9:01, Jérôme Hénin 写道:
>
> Dear Haochuan,
>
> For your purpose the simplest approach would seem to define multiple RMSD
> components and process them in a script, similar to the patch collective
> variables.
>
> You could do something similar in C++, but then I would use a vector of
> RMSD components, rather than a vector of atom groups.
>
> Best,
> Jerome
>
> On Sat, 9 Mar 2019 at 02:44, yjcoshc <yjcoshc_at_gmail.com> wrote:
>
>> Dear NAMD and Colvars developers,
>>
>> I am trying to implement a collective variable (CV) through the C++
>> interface of Colvars. This CV should accept multiple reference frames
>> and compute the minimal RMSDs of the current atomic coordinates with
>> respect of them. In the class of new CV I have tried to add a vector
>> like "std::vector<cvm::atom_group*> comp_atoms" and then parsed the
>> selected atoms in the constructor like this loop:
>>
>> for (size_t i_frame = 0; i_frame < reference_frames.size();
>> ++i_frame) {
>> cvm::atom_group* tmp_atoms = parse_group(conf, "atoms");
>> // Swipe from the rmsd class
>> tmp_atoms->b_center = true;
>> tmp_atoms->b_rotate = true;
>> tmp_atoms->ref_pos = reference_frames[i_frame];
>> tmp_atoms->center_ref_pos();
>> tmp_atoms->disable(f_ag_fit_gradients);
>> tmp_atoms->rot.request_group1_gradients(tmp_atoms->size());
>> tmp_atoms->rot.request_group2_gradients(tmp_atoms->size());
>> comp_atoms.push_back(tmp_atoms);
>> }
>>
>> However this doesn't work, and the colvars ends up with throwing an
>> exception of "std::out_of_range" in "basic_string::erase" after the
>> constructor finishes.
>>
>> Is there a standard way to compute multiple minimal RMSDs in the C++
>> interface?
>>
>> Thanks,
>>
>> Haochuan Chen
>>
>>

-- 
Giacomo Fiorin
Associate Professor of Research, Temple University, Philadelphia, PA
Contractor, National Institutes of Health, Bethesda, MD
http://goo.gl/Q3TBQU
https://github.com/giacomofiorin

This archive was generated by hypermail 2.1.6 : Tue Dec 31 2019 - 23:20:35 CST