Skip to content

Commit

Permalink
User-defined index group names, add dihedralPC
Browse files Browse the repository at this point in the history
  • Loading branch information
@jhenin
jhenin committed Apr 5, 2024
1 parent 01c70ed commit cb738d8
Showing 1 changed file with 72 additions and 37 deletions.
109 changes: 72 additions & 37 deletions src/colvarcomp_protein.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ int colvar::alpha_angles::init(std::string const &conf)
cvm::atom_group group_CA, group_N, group_O;

std::string residues_conf = "";
std::string prefix;

// residueRange is mandatory for the topology-based case
if (key_lookup(conf, "residueRange", &residues_conf)) {
Expand All @@ -62,19 +63,21 @@ int colvar::alpha_angles::init(std::string const &conf)

} else {
b_use_index_groups = true;
get_keyval(conf, "prefix", prefix, "alpha");

// Not all groups are mandatory, parse silently
group_CA.add_index_group("alphaCA", true);
group_N.add_index_group("alphaN", true);
group_O.add_index_group("alphaO", true);
group_CA.add_index_group(prefix + "CA", true);
group_N.add_index_group(prefix + "N", true);
group_O.add_index_group(prefix + "O", true);
int na = group_CA.size();
int nn = group_N.size();
int no = group_O.size();
if ((nn != 0 || no != 0) && (nn != no)) {
return cvm::error("Error: If either is provided, atom groups alphaN and alphaO must have the same number of atoms.",
return cvm::error("Error: If either is provided, atom groups " + prefix + "N and " + prefix + "O must have the same number of atoms.",
COLVARS_INPUT_ERROR);
}
if (nn != 0 && na != 0 && nn != na) {
return cvm::error("Error: If both are provided, atom groups alphaN and alphaCA must have the same number of atoms.",
return cvm::error("Error: If both are provided, atom groups " + prefix + "N and " + prefix + "CA must have the same number of atoms.",
COLVARS_INPUT_ERROR);
}
}
Expand All @@ -95,7 +98,7 @@ int colvar::alpha_angles::init(std::string const &conf)
if (hb_coeff < 1.0) {
if (b_use_index_groups) {
if (group_CA.size() < 5) {
return cvm::error("Not enough atoms (" + cvm::to_str(group_CA.size()) + ") in index group \"alphaCA\"",
return cvm::error("Not enough atoms (" + cvm::to_str(group_CA.size()) + ") in index group \"" + prefix + "CA\"",
COLVARS_INPUT_ERROR);
}
for (size_t i = 0; i < group_CA.size()-2; i++) {
Expand Down Expand Up @@ -130,7 +133,7 @@ int colvar::alpha_angles::init(std::string const &conf)
if (hb_coeff > 0.0) {
if (b_use_index_groups) {
if (group_N.size() < 5) {
return cvm::error("Not enough atoms (" + cvm::to_str(group_N.size()) + ") in index group \"alphaN\"",
return cvm::error("Not enough atoms (" + cvm::to_str(group_N.size()) + ") in index group \"" + prefix + "N\"",
COLVARS_INPUT_ERROR);
}
for (size_t i = 0; i < group_N.size()-4; i++) {
Expand Down Expand Up @@ -338,41 +341,62 @@ int colvar::dihedPC::init(std::string const &conf)
if (cvm::debug())
cvm::log("Initializing dihedral PC object.\n");

bool b_use_index_groups = false;
std::string segment_id;
get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));

std::vector<int> residues;
{
std::string residues_conf = "";
key_lookup(conf, "residueRange", &residues_conf);
size_t n_residues;
std::string residues_conf = "";
std::string prefix;
cvm::atom_group group_CA, group_N, group_C;

// residueRange is mandatory for the topology-based case
if (key_lookup(conf, "residueRange", &residues_conf)) {
if (residues_conf.size()) {
std::istringstream is(residues_conf);
int initial, final;
char dash;
if ( (is >> initial) && (initial > 0) &&
(is >> dash) && (dash == '-') &&
(is >> final) && (final > 0) ) {
(is >> dash) && (dash == '-') &&
(is >> final) && (final > 0) ) {
for (int rnum = initial; rnum <= final; rnum++) {
residues.push_back(rnum);
}
}
} else {
error_code |=
cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
return cvm::error("Error: no residues defined in \"residueRange\".\n", COLVARS_INPUT_ERROR);
}
}
n_residues = residues.size();
get_keyval(conf, "psfSegID", segment_id, std::string("MAIN"));

} else {

if (residues.size() < 2) {
b_use_index_groups = true;
get_keyval(conf, "prefix", prefix, "dihed");

// Not all groups are mandatory, parse silently
group_CA.add_index_group(prefix + "CA", true);
group_N.add_index_group(prefix + "N", true);
group_C.add_index_group(prefix + "C", true);
int na = group_CA.size();
int nn = group_N.size();
int nc = group_C.size();
if ((nn != na || na != nc)) {
return cvm::error("Error: atom groups " + prefix + "N, " + prefix + "CA, and " + prefix +
"C must have the same number of atoms.", COLVARS_INPUT_ERROR);
}
n_residues = nn;
}
if (n_residues < 2) {
error_code |=
cvm::error("Error: dihedralPC requires at least two residues.\n", COLVARS_INPUT_ERROR);
cvm::error("Error: dihedralPC requires at least two residues.\n", COLVARS_INPUT_ERROR);
}

std::string const &sid = segment_id;
std::vector<int> const &r = residues;

std::string vecFileName;
int vecNumber;
if (get_keyval(conf, "vectorFile", vecFileName, vecFileName)) {
int vecNumber;
get_keyval(conf, "vectorNumber", vecNumber, 0);
if (vecNumber < 1) {
error_code |=
Expand All @@ -387,9 +411,8 @@ int colvar::dihedPC::init(std::string const &conf)
}

// TODO: adapt to different formats by setting this flag
bool eigenvectors_as_columns = true;

if (eigenvectors_as_columns) {
// bool eigenvectors_as_columns = true;
// if (eigenvectors_as_columns) {
// Carma-style dPCA file
std::string line;
cvm::real c;
Expand All @@ -400,9 +423,7 @@ int colvar::dihedPC::init(std::string const &conf)
for (int i=0; i<vecNumber; i++) ls >> c;
coeffs.push_back(c);
}
}
/* TODO Uncomment this when different formats are recognized
else {
/* } else { // Uncomment this when different formats are recognized
// Eigenvectors as lines
// Skip to the right line
for (int i = 1; i<vecNumber; i++)
Expand All @@ -428,28 +449,42 @@ int colvar::dihedPC::init(std::string const &conf)
get_keyval(conf, "vector", coeffs, coeffs);
}

if ( coeffs.size() != 4 * (residues.size() - 1)) {
if ( coeffs.size() != 4 * (n_residues - 1)) {
error_code |= cvm::error("Error: wrong number of coefficients: " + cvm::to_str(coeffs.size()) +
". Expected " + cvm::to_str(4 * (residues.size() - 1)) +
". Expected " + cvm::to_str(4 * (n_residues - 1)) +
" (4 coeffs per residue, minus one residue).\n",
COLVARS_INPUT_ERROR);
}

for (size_t i = 0; i < residues.size()-1; i++) {
for (size_t i = 0; i < n_residues-1; i++) {
// Psi
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "N", sid),
cvm::atom(r[i ], "CA", sid),
cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid)));
if (b_use_index_groups) {
theta.push_back(new colvar::dihedral( group_N[i],
group_CA[i],
group_C[i],
group_N[i+1]));
} else {
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "N", sid),
cvm::atom(r[i ], "CA", sid),
cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid)));
}
register_atom_group(theta.back()->atom_groups[0]);
register_atom_group(theta.back()->atom_groups[1]);
register_atom_group(theta.back()->atom_groups[2]);
register_atom_group(theta.back()->atom_groups[3]);
// Phi (next res)
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid),
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+1], "C", sid)));
if (b_use_index_groups) {
theta.push_back(new colvar::dihedral(group_C[i],
group_N[i+1],
group_CA[i+1],
group_C[i+1]));
} else {
theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "C", sid),
cvm::atom(r[i+1], "N", sid),
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+1], "C", sid)));
}
register_atom_group(theta.back()->atom_groups[0]);
register_atom_group(theta.back()->atom_groups[1]);
register_atom_group(theta.back()->atom_groups[2]);
Expand Down

0 comments on commit cb738d8

Please sign in to comment.