Homec4science
Diffusion lammps (lammm-master) Browse

lammps/lib/colvars/colvarcomp_protein.cpp
fee52f2fac2blammm-master

Last Change
Raw File

colvarcomp_protein.cpp
View Options

#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
//////////////////////////////////////////////////////////////////////
// alpha component
//////////////////////////////////////////////////////////////////////
colvar::alpha_angles::alpha_angles (std::string const &conf)
: cvc (conf)
{
if (cvm::debug())
cvm::log ("Initializing alpha_angles object.\n");
function_type = "alpha_angles";
x.type (colvarvalue::type_scalar);
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);
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) ) {
for (int rnum = initial; rnum <= final; rnum++) {
residues.push_back (rnum);
}
}
} else {
cvm::fatal_error ("Error: no residues defined in \"residueRange\".\n");
}
}
if (residues.size() < 5) {
cvm::fatal_error ("Error: not enough residues defined in \"residueRange\".\n");
}
std::string const &sid = segment_id;
std::vector<int> const &r = residues;
get_keyval (conf, "hBondCoeff", hb_coeff, 0.5);
if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) {
cvm::fatal_error ("Error: hBondCoeff must be defined between 0 and 1.\n");
}
get_keyval (conf, "angleRef", theta_ref, 88.0);
get_keyval (conf, "angleTol", theta_tol, 15.0);
if (hb_coeff < 1.0) {
for (size_t i = 0; i < residues.size()-2; i++) {
theta.push_back (new colvar::angle (cvm::atom (r[i ], "CA", sid),
cvm::atom (r[i+1], "CA", sid),
cvm::atom (r[i+2], "CA", sid)));
}
} else {
cvm::log ("The hBondCoeff specified will disable the Calpha-Calpha-Calpha angle terms.\n");
}
{
cvm::real r0;
size_t en, ed;
get_keyval (conf, "hBondCutoff", r0, (3.3 * cvm::unit_angstrom()));
get_keyval (conf, "hBondExpNumer", en, 6);
get_keyval (conf, "hBondExpDenom", ed, 8);
if (hb_coeff > 0.0) {
for (size_t i = 0; i < residues.size()-4; i++) {
hb.push_back (new colvar::h_bond (cvm::atom (r[i ], "O", sid),
cvm::atom (r[i+4], "N", sid),
r0, en, ed));
}
} else {
cvm::log ("The hBondCoeff specified will disable the hydrogen bond terms.\n");
}
}
if (cvm::debug())
cvm::log ("Done initializing alpha_angles object.\n");
}
colvar::alpha_angles::alpha_angles()
: cvc ()
{
function_type = "alpha_angles";
x.type (colvarvalue::type_scalar);
}
colvar::alpha_angles::~alpha_angles()
{
while (theta.size() != 0) {
delete theta.back();
theta.pop_back();
}
while (hb.size() != 0) {
delete hb.back();
hb.pop_back();
}
}
void colvar::alpha_angles::calc_value()
{
x.real_value = 0.0;
if (theta.size()) {
cvm::real const theta_norm =
(1.0-hb_coeff) / cvm::real (theta.size());
for (size_t i = 0; i < theta.size(); i++) {
(theta[i])->calc_value();
cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
cvm::real const f = ( (1.0 - std::pow (t, (int) 2)) /
(1.0 - std::pow (t, (int) 4)) );
x.real_value += theta_norm * f;
if (cvm::debug())
cvm::log ("Calpha-Calpha angle no. "+cvm::to_str (i+1)+" in \""+
this->name+"\" has a value of "+
(cvm::to_str ((theta[i])->value().real_value))+
" degrees, f = "+cvm::to_str (f)+".\n");
}
}
if (hb.size()) {
cvm::real const hb_norm =
hb_coeff / cvm::real (hb.size());
for (size_t i = 0; i < hb.size(); i++) {
(hb[i])->calc_value();
x.real_value += hb_norm * (hb[i])->value().real_value;
if (cvm::debug())
cvm::log ("Hydrogen bond no. "+cvm::to_str (i+1)+" in \""+
this->name+"\" has a value of "+
(cvm::to_str ((hb[i])->value().real_value))+".\n");
}
}
}
void colvar::alpha_angles::calc_gradients()
{
for (size_t i = 0; i < theta.size(); i++)
(theta[i])->calc_gradients();
for (size_t i = 0; i < hb.size(); i++)
(hb[i])->calc_gradients();
}
void colvar::alpha_angles::apply_force (colvarvalue const &force)
{
if (theta.size()) {
cvm::real const theta_norm =
(1.0-hb_coeff) / cvm::real (theta.size());
for (size_t i = 0; i < theta.size(); i++) {
cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
cvm::real const f = ( (1.0 - std::pow (t, (int) 2)) /
(1.0 - std::pow (t, (int) 4)) );
cvm::real const dfdt =
1.0/(1.0 - std::pow (t, (int) 4)) *
( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow (t, (int) 3)) );
(theta[i])->apply_force (theta_norm *
dfdt * (1.0/theta_tol) *
force.real_value );
}
}
if (hb.size()) {
cvm::real const hb_norm =
hb_coeff / cvm::real (hb.size());
for (size_t i = 0; i < hb.size(); i++) {
(hb[i])->apply_force (0.5 * hb_norm * force.real_value);
}
}
}
//////////////////////////////////////////////////////////////////////
// dihedral principal component
//////////////////////////////////////////////////////////////////////
// FIXME: this will not make collect_gradients work
// because gradients in individual atom groups
// are those of the sub-cvcs (angle, hb), not those
// of this cvc (alpha)
// This is true of all cvcs with sub-cvcs, and those
// that do not calculate explicit gradients
// SO: we need a flag giving the availability of
// atomic gradients
colvar::dihedPC::dihedPC (std::string const &conf)
: cvc (conf)
{
if (cvm::debug())
cvm::log ("Initializing dihedral PC object.\n");
function_type = "dihedPC";
x.type (colvarvalue::type_scalar);
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);
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) ) {
for (int rnum = initial; rnum <= final; rnum++) {
residues.push_back (rnum);
}
}
} else {
cvm::fatal_error ("Error: no residues defined in \"residueRange\".\n");
}
}
if (residues.size() < 2) {
cvm::fatal_error ("Error: dihedralPC requires at least two residues.\n");
}
std::string const &sid = segment_id;
std::vector<int> const &r = residues;
std::string vecFileName;
int vecNumber;
if (get_keyval (conf, "vectorFile", vecFileName, vecFileName)) {
get_keyval (conf, "vectorNumber", vecNumber, 0);
if (vecNumber < 1)
cvm::fatal_error ("A positive value of vectorNumber is required.");
std::ifstream vecFile;
vecFile.open (vecFileName.c_str());
if (!vecFile.good())
cvm::fatal_error ("Error opening dihedral PCA vector file " + vecFileName + " for reading");
// TODO: adapt to different formats by setting this flag
bool eigenvectors_as_columns = true;
if (eigenvectors_as_columns) {
// Carma-style dPCA file
std::string line;
cvm::real c;
while (vecFile.good()) {
getline(vecFile, line);
if (line.length() < 2) break;
std::istringstream ls (line);
for (int i=0; i<vecNumber; i++) ls >> c;
coeffs.push_back(c);
}
}
/* TODO Uncomment this when different formats are recognized
else {
// Eigenvectors as lines
// Skip to the right line
for (int i = 1; i<vecNumber; i++)
vecFile.ignore(999999, '\n');
if (!vecFile.good())
cvm::fatal_error ("Error reading dihedral PCA vector file " + vecFileName);
std::string line;
getline(vecFile, line);
std::istringstream ls (line);
cvm::real c;
while (ls.good()) {
ls >> c;
coeffs.push_back(c);
}
}
*/
vecFile.close();
} else {
get_keyval (conf, "vector", coeffs, coeffs);
}
if ( coeffs.size() != 4 * (residues.size() - 1)) {
cvm::fatal_error ("Error: wrong number of coefficients: " +
cvm::to_str(coeffs.size()) + ". Expected " +
cvm::to_str(4 * (residues.size() - 1)) +
" (4 coeffs per residue, minus one residue).\n");
}
for (size_t i = 0; i < residues.size()-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)));
// 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 (cvm::debug())
cvm::log ("Done initializing dihedPC object.\n");
}
colvar::dihedPC::dihedPC()
: cvc ()
{
function_type = "dihedPC";
x.type (colvarvalue::type_scalar);
}
colvar::dihedPC::~dihedPC()
{
while (theta.size() != 0) {
delete theta.back();
theta.pop_back();
}
}
void colvar::dihedPC::calc_value()
{
x.real_value = 0.0;
for (size_t i = 0; i < theta.size(); i++) {
theta[i]->calc_value();
cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
x.real_value += coeffs[2*i ] * std::cos(t)
+ coeffs[2*i+1] * std::sin(t);
}
}
void colvar::dihedPC::calc_gradients()
{
for (size_t i = 0; i < theta.size(); i++) {
theta[i]->calc_gradients();
}
}
void colvar::dihedPC::apply_force (colvarvalue const &force)
{
for (size_t i = 0; i < theta.size(); i++) {
cvm::real const t = (PI / 180.) * theta[i]->value().real_value;
cvm::real const dcosdt = - (PI / 180.) * std::sin(t);
cvm::real const dsindt = (PI / 180.) * std::cos(t);
theta[i]->apply_force((coeffs[2*i ] * dcosdt +
coeffs[2*i+1] * dsindt) * force);
}
}
c4science · Help