Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[RFC] Work towards reusable CVCs #700

Draft
wants to merge 5 commits into
base: reusable-cvcs
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 149 additions & 0 deletions namd/src/colvarproxy_namd.C
Original file line number Diff line number Diff line change
Expand Up @@ -41,10 +41,19 @@
#include "colvarproxy.h"
#include "colvarproxy_namd.h"
#include "colvarscript.h"
#include "colvarcomp.h"

#ifdef USE_CKLOOP
#include "CkLambda.h"
#endif



colvarproxy_namd::colvarproxy_namd()
{
#if defined(CMK_SMP)
charm_lock_state = CmiCreateLock();
#endif
engine_name_ = "NAMD";

version_int = get_version_from_string(COLVARPROXY_VERSION);
Expand Down Expand Up @@ -151,6 +160,9 @@ colvarproxy_namd::colvarproxy_namd()

colvarproxy_namd::~colvarproxy_namd()
{
#if defined(CMK_SMP)
CmiDestroyLock(charm_lock_state);
#endif
delete reduction;
}

Expand Down Expand Up @@ -1539,6 +1551,143 @@ int colvarproxy_namd::smp_colvars_loop()
return cvm::get_error();
}

struct cvc_info {
bool enabled;
bool parent_colvar_total_force_calc;
bool parent_colvar_Jacobian;
std::shared_ptr<colvar::cvc> cvc_ptr;
};

// TODO: Too bad! We have to traverse the AST twice!
void cvc_max_depth(
const std::vector<std::shared_ptr<colvar::cvc>>& cvcs,
std::unordered_map<std::shared_ptr<colvar::cvc>, int>& cvc_depth_map,
int current_level = 0) {
for (auto it = cvcs.begin(); it != cvcs.end(); ++it) {
auto map_it = cvc_depth_map.find(*it);
if (map_it == cvc_depth_map.end()) {
cvc_depth_map.insert({*it, current_level});
} else {
if (map_it->second < current_level) {
map_it->second = current_level;
}
}
cvc_max_depth((*it)->children_cvcs(), cvc_depth_map, current_level+1);
}
}

void flatten_all_cvc_tree(
std::shared_ptr<colvar::cvc>& parent,
std::unordered_map<std::shared_ptr<colvar::cvc>, cvc_info>& cvc_info_map,
bool parent_colvar_total_force_calc,
bool parent_colvar_Jacobian) {
auto children = parent->children_cvcs();
for (auto it = children.begin(); it != children.end(); ++it) {
flatten_all_cvc_tree(*it, cvc_info_map, parent_colvar_total_force_calc, parent_colvar_Jacobian);
}
auto it_find = cvc_info_map.find(parent);
if (it_find == cvc_info_map.end()) {
cvc_info_map.insert({parent, cvc_info{
// TODO: Here calc_cvc_values calls cvcs[i]->is_enabled() , which is is_enabled(int f = f_cv_active)
// I know both f_cv_active and f_cvc_active are 0 but are they the same option??
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The fact that f_cv_active and f_cvc_active have the same numerical value has no consequence, they should never be used in the same context, because they are respectively only meaningful in a colvar or cvc object. The colvardeps data of these two classes are non-overlapping. The relationship between them is a vertical (parent-child) dependency.

If we were to merge the two levels, then these two features would merge.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The issue is that in

if (!cvcs[i]->is_enabled()) continue;
, is_enabled() is called, and since there is no function parameter passed, so I think it would call
inline bool is_enabled(int f = f_cv_active) const {

which checks f_cv_active instead of f_cvc_active. My code is to follow what the original calls in calc_cvc_values and I am not sure if I need to follow it the same way.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You're right, sorry. That is somewhat sloppy writing that I didn't remember well. It does rely on the "active" property being number 0.

// NOTE that all feature enums should start with f_*_active

But a class inheriting from colvardeps could also override is_enabled() and change this convention.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How much effort do you think it would be to remove the default argument for this virtual function? Remember that this is one of the "issues" that clang-tidy was complaining about.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That would be very little effort. I'm happy to do that if that helps in any way.

parent->is_enabled(colvardeps::features_cvc::f_cvc_active),
parent_colvar_total_force_calc, parent_colvar_Jacobian,
parent
}});
}
}

int colvarproxy_namd::smp_colvars_loop2() {
colvarmodule *cv = this->colvars;
// Bypass the colvar objects and find all CVC objects
std::vector<std::shared_ptr<colvar::cvc>> all_cvcs;
std::unordered_map<std::shared_ptr<colvar::cvc>, cvc_info> cvc_info_map;
for (auto it = cv->variables_active()->begin();
it != cv->variables_active()->end(); ++it) {
// TODO: Bad design! What will happen if CVC a is in a "colvar" block
// that does not support total_force_calc, but is then reused in
// another block that requires total_force_calc even if it supports Jacobian itself???
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see the problem here. Assuming a cvc can have several parents: the colvar that does require a total force can enable it in its children cvcs, even if other parents don't require (or support) it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right. I just thought that if a colvar does not require the total force, then it would disable the corresponding feature of the children, but it seems the code do not check the dependency in that way.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's the other way: disabled by default, and enabled on request - then disabled again if the refcount falls to zero.

const bool total_force_on = (*it)->is_enabled(colvardeps::features_colvar::f_cv_total_force_calc);
const bool Jacobian_on = (*it)->is_enabled(colvardeps::features_colvar::f_cv_Jacobian);
// use_total_forces.push_back(total_force_on);
for (auto it_cvc = (*it)->cvcs_begin();
it_cvc != (*it)->cvcs_end(); ++it_cvc) {
std::shared_ptr<colvar::cvc> p = *it_cvc;
flatten_all_cvc_tree(p, cvc_info_map, total_force_on, Jacobian_on);
all_cvcs.push_back(p);
}
}
// Since some CVCs may depend on others, the CVCs constitutes a directed acyclic graph.
// Walk through the graph, and determine the maximum depth of each CVC.
std::unordered_map<std::shared_ptr<colvar::cvc>, int> cvc_depth_map;
cvc_max_depth(all_cvcs, cvc_depth_map, 0);
// Determine the max depth of the graph
int max_depth = 0;
for (auto it = cvc_depth_map.begin(); it != cvc_depth_map.end(); ++it) {
if (it->second > max_depth) max_depth = it->second;
}
// Group the CVCs according to their max depths
// TODO: Can we save this vector since the CVCs are less likely changed over time?
// But what should I do if the user changes the Colvars config dynamically? I wish there is a signal/slot...
std::vector<std::vector<cvc_info>> sorted_cvcs(max_depth+1);
for (auto it = cvc_depth_map.begin(); it != cvc_depth_map.end(); ++it) {
sorted_cvcs[it->second].push_back(cvc_info_map[it->first]);
}
const bool tf_switch_1 = (cvm::step_relative() > 0) && (!total_forces_same_step());
const bool tf_switch_2 = total_forces_same_step();
// Calculate the CVCs in parallel. The high-depth ones should be calculated at first.
cvm::increase_depth();
for (auto it = sorted_cvcs.rbegin(); it != sorted_cvcs.rend(); ++it) {
const int numChunks = it->size();
const int lowerRange = 0;
const int upperRange = numChunks - 1;
cvc_info* data = it->data();
auto lambda_fn = [data, tf_switch_1, tf_switch_2](int start, int end, void* unused){
for (int i = start; i <= end; i++) {
std::shared_ptr<colvar::cvc>& cvc_ptr = data[i].cvc_ptr;
if (data[i].enabled) {
// Follow the order of colvar::calc_cvcs:
// 1. calc_cvc_total_force
if (data[i].parent_colvar_total_force_calc && tf_switch_1) {
cvc_ptr->calc_force_invgrads();
}
// 2. calc_cvc_values
// TODO: Oops! read_data() clears the gradients!
// CVC a has a sub-CVC b, and we run b->calc_gradients at
// first and followed by a->read_data, then the previous gradients are cleared!!
// Without major refactoring, it seems the only way to avoid the problem is
// not to register the atom groups of sub-CVCs to the parent CVC.
// Then I have to rely on modify_children_cvcs_atom_gradients to do backward propagation.
// PLUMED has PLMD::ActionAtomistic and PLMD::ActionWithValue but Colvars has none.
cvc_ptr->read_data();
cvc_ptr->calc_value();
// 3. calc_cvc_gradients
if (cvc_ptr->is_enabled(colvardeps::features_cvc::f_cvc_gradient)) {
cvc_ptr->calc_gradients();
cvc_ptr->calc_fit_gradients();
// TODO: Implement debug_gradients for CVC of sub-CVCs
if (cvc_ptr->is_enabled(colvardeps::features_cvc::f_cvc_debug_gradient)) {
cvc_ptr->debug_gradients();
}
}
// 4. calc_cvc_Jacobians
if (data[i].parent_colvar_Jacobian) {
cvc_ptr->calc_Jacobian_derivative();
}
// 5. calc_cvc_total_force
if (data[i].parent_colvar_total_force_calc && tf_switch_2) {
cvc_ptr->calc_force_invgrads();
}
}
}
};
CkLoop_Parallelize(numChunks, lowerRange, upperRange,
lambda_fn, NULL, CKLOOP_NONE, NULL);
}
cvm::decrease_depth();
return cvm::get_error();
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Parallelization and checking the error

When I wrote and test new code, I found that NAMD may not exit cleanly (some other threads were still running) if one thread was exited. I think the implementation of colvarproxy_namd::error or at least its use in the parallel region, needs to be revised to use a std::atomic<bool> to setup an error "bit" instead of calling NAMD_die directly. Then after all threads are finished and joined, we can use cvm::get_error() to check the error bit and exit if there is an error.

}


void calc_cv_biases_smp(int first, int last, void *result, int paramNum, void *param)
{
Expand Down
14 changes: 11 additions & 3 deletions namd/src/colvarproxy_namd.h
Original file line number Diff line number Diff line change
Expand Up @@ -123,6 +123,7 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster {
int check_smp_enabled() override;

int smp_colvars_loop() override;
int smp_colvars_loop2() override;

int smp_biases_loop();

Expand Down Expand Up @@ -150,18 +151,25 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster {

int smp_lock()
{
charm_lock_state = CmiCreateLock();
CmiLock(charm_lock_state);
// TODO: Why does the old code create lock instead of holding the existing lock???
// charm_lock_state = CmiCreateLock();
return COLVARS_OK;
}

int smp_trylock()
{
return COLVARS_NOT_IMPLEMENTED;
const int ret = CmiTryLock(charm_lock_state);
if (ret == 0) return COLVARS_OK;
else return COLVARS_ERROR;
// return COLVARS_NOT_IMPLEMENTED;
}

int smp_unlock()
{
CmiDestroyLock(charm_lock_state);
// NOTE: It seems the smp locking was never implemented correctly. No matter I got so many mysterious crashes with it!
// CmiDestroyLock(charm_lock_state);
CmiUnlock(charm_lock_state);
return COLVARS_OK;
}

Expand Down
Loading