Bugfix for Colvars library (update to version 2019-08-05)

Bugfix for group2CenterOnly (coordNum option):
https://github.com/Colvars/colvars/pull/278
This commit is contained in:
Giacomo Fiorin
2019-08-05 14:12:38 -04:00
parent 096c225594
commit 51ba9bd520
5 changed files with 80 additions and 97 deletions

View File

@ -91,7 +91,7 @@ cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
colvar::coordnum::coordnum(std::string const &conf)
: cvc(conf), b_anisotropic(false), group2_center(NULL), pairlist(NULL)
: cvc(conf), b_anisotropic(false), pairlist(NULL)
{
function_type = "coordnum";
@ -156,14 +156,7 @@ colvar::coordnum::coordnum(std::string const &conf)
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
bool b_group2_center_only = false;
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
if (b_group2_center_only) {
if (!group2_center) {
group2_center = new cvm::atom_group();
group2_center->add_atom(cvm::atom());
}
}
get_keyval(conf, "tolerance", tolerance, 0.0);
if (tolerance > 0) {
@ -185,7 +178,7 @@ colvar::coordnum::coordnum(std::string const &conf)
colvar::coordnum::coordnum()
: b_anisotropic(false), group2_center(NULL), pairlist(NULL)
: b_anisotropic(false), pairlist(NULL)
{
function_type = "coordnum";
x.type(colvarvalue::type_scalar);
@ -197,69 +190,60 @@ colvar::coordnum::~coordnum()
if (pairlist != NULL) {
delete [] pairlist;
}
if (group2_center != NULL) {
delete group2_center;
}
template<int flags> void colvar::coordnum::main_loop(bool **pairlist_elem)
{
if (b_group2_center_only) {
cvm::atom group2_com_atom;
group2_com_atom.pos = group2->center_of_mass();
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, group2_com_atom,
pairlist_elem,
tolerance);
}
if (b_group2_center_only) {
group2->set_weighted_gradient(group2_com_atom.grad);
}
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++) {
for (cvm::atom_iter ai2 = group2->begin(); ai2 != group2->end(); ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
pairlist_elem,
tolerance);
}
}
}
}
template<int compute_flags> int colvar::coordnum::compute_coordnum()
{
if (group2_center) {
(*group2_center)[0].pos = group2->center_of_mass();
group2_center->calc_required_properties();
}
cvm::atom_group *group2p = group2_center ? group2_center : group2;
bool const use_pairlist = (pairlist != NULL);
bool const rebuild_pairlist = (pairlist != NULL) &&
(cvm::step_relative() % pairlist_freq == 0);
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
cvm::atom_iter ai1 = group1->begin(), ai2 = group2p->begin();
cvm::atom_iter const ai1_end = group1->end();
cvm::atom_iter const ai2_end = group2p->end();
if (b_anisotropic) {
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
main_loop<flags>(&pairlist_elem);
} else {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
main_loop<flags>(&pairlist_elem);
}
} else { // if (use_pairlist) {
} else {
int const flags = compute_flags | ef_anisotropic;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
main_loop<flags>(NULL);
}
} else {
@ -267,46 +251,17 @@ template<int compute_flags> int colvar::coordnum::compute_coordnum()
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
main_loop<flags>(&pairlist_elem);
} else {
int const flags = compute_flags | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
main_loop<flags>(&pairlist_elem);
}
} else { // if (use_pairlist) {
} else {
int const flags = compute_flags;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
}
}
if (compute_flags & ef_gradients) {
if (group2_center) {
group2->set_weighted_gradient((*group2_center)[0].grad);
main_loop<flags>(NULL);
}
}