ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

update colvars module with support for index files

Browse Source
This commit is contained in:
Axel Kohlmeyer
2013-05-15 15:17:47 +02:00
parent 92b0626589
commit f43a004f01
7 changed files with 261 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

19
lib/colvars/colvaratoms.cpp
View File

@ -129,6 +129,25 @@ void cvm::atom_group::parse (std::string const &conf,
atom_indexes.clear();
}
std::string index_group_name;
if (get_keyval (group_conf, "indexGroup", index_group_name)) {
// use an index group from the index file read globally
std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
for ( ; names_i != cvm::index_group_names.end() ; names_i++, index_groups_i++) {
if (*names_i == index_group_name)
break;
}
if (names_i == cvm::index_group_names.end()) {
cvm::fatal_error ("Error: could not find index group "+
index_group_name+" among those provided by the index file.\n");
}
this->reserve (index_groups_i->size());
for (size_t i = 0; i < index_groups_i->size(); i++) {
this->push_back (cvm::atom ((*index_groups_i)[i]));
}
}
}
{

122
lib/colvars/colvarbias.cpp
View File

@ -41,10 +41,11 @@ colvarbias::colvarbias (std::string const &conf, char const *key)
add_colvar (colvars_str[i]);
}
}
if (!colvars.size()) {
cvm::fatal_error ("Error: no collective variables specified.\n");
}
get_keyval (conf, "outputEnergy", b_output_energy, false);
}
@ -82,6 +83,40 @@ void colvarbias::communicate_forces()
}
void colvarbias::change_configuration(std::string const &conf)
{
cvm::fatal_error ("Error: change_configuration() not implemented.\n");
}
cvm::real colvarbias::energy_difference(std::string const &conf)
{
cvm::fatal_error ("Error: energy_difference() not implemented.\n");
return 0.;
}
std::ostream & colvarbias::write_traj_label (std::ostream &os)
{
os << " ";
if (b_output_energy)
os << " E_"
<< cvm::wrap_string (this->name, cvm::en_width-2);
return os;
}
std::ostream & colvarbias::write_traj (std::ostream &os)
{
os << " ";
if (b_output_energy)
os << " "
<< bias_energy;
return os;
}
colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
char const *key)
@ -171,25 +206,16 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
}
}
get_keyval (conf, "outputCenters", b_output_centers, false);
get_keyval (conf, "outputAccumulatedWork", b_output_acc_work, false);
acc_work = 0.0;
if (cvm::debug())
cvm::log ("Done initializing a new harmonic restraint bias.\n");
}
void colvarbias::change_configuration(std::string const &conf)
{
cvm::fatal_error ("Error: change_configuration() not implemented.\n");
}
cvm::real colvarbias::energy_difference(std::string const &conf)
{
cvm::fatal_error ("Error: energy_difference() not implemented.\n");
return 0.;
}
void colvarbias_harmonic::change_configuration(std::string const &conf)
void colvarbias_harmonic::change_configuration (std::string const &conf)
{
get_keyval (conf, "forceConstant", force_k, force_k);
if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
@ -379,6 +405,15 @@ cvm::real colvarbias_harmonic::update()
colvar_centers[i]))+"\n");
}
if (b_output_acc_work) {
if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) {
for (size_t i = 0; i < colvars.size(); i++) {
// project forces on the calculated increments at this step
acc_work += colvar_forces[i] * centers_incr[i];
}
}
}
if (cvm::debug())
cvm::log ("Current forces for the harmonic bias \""+
this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
@ -442,6 +477,11 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is)
cvm::fatal_error ("Error: current stage is missing from the restart.\n");
}
if (b_output_acc_work) {
if (!get_keyval (conf, "accumulatedWork", acc_work))
cvm::fatal_error ("Error: accumulatedWork is missing from the restart.\n");
}
is >> brace;
if (brace != "}") {
cvm::fatal_error ("Error: corrupt restart information for harmonic bias \""+
@ -484,9 +524,61 @@ std::ostream & colvarbias_harmonic::write_restart (std::ostream &os)
<< stage << "\n";
}
if (b_output_acc_work) {
os << " accumulatedWork " << acc_work << "\n";
}
os << " }\n"
<< "}\n\n";
return os;
}
std::ostream & colvarbias_harmonic::write_traj_label (std::ostream &os)
{
os << " ";
if (b_output_energy)
os << " E_"
<< cvm::wrap_string (this->name, cvm::en_width-2);
if (b_output_centers)
for (size_t i = 0; i < colvars.size(); i++) {
size_t const this_cv_width = (colvars[i]->value()).output_width (cvm::cv_width);
os << " x0_"
<< cvm::wrap_string (colvars[i]->name, this_cv_width-3);
}
if (b_output_acc_work)
os << " W_"
<< cvm::wrap_string (this->name, cvm::en_width-2);
return os;
}
std::ostream & colvarbias_harmonic::write_traj (std::ostream &os)
{
os << " ";
if (b_output_energy)
os << " "
<< std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
<< bias_energy;
if (b_output_centers)
for (size_t i = 0; i < colvars.size(); i++) {
os << " "
<< std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
<< colvar_centers[i];
}
if (b_output_acc_work)
os << " "
<< std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
<< acc_work;
return os;
}

80
lib/colvars/colvarbias.h
View File

@ -52,6 +52,13 @@ public:
/// Write the bias configuration to a restart file
virtual std::ostream & write_restart (std::ostream &os) = 0;
/// Write a label to the trajectory file (comment line)
virtual std::ostream & write_traj_label (std::ostream &os);
/// Output quantities such as the bias energy to the trajectory file
virtual std::ostream & write_traj (std::ostream &os);
protected:
/// \brief Pointers to collective variables to which the bias is
@ -62,10 +69,12 @@ protected:
/// \brief Current forces from this bias to the colvars
std::vector<colvarvalue> colvar_forces;
/// \brief Current energy of this bias (colvar_forces should be
/// obtained by deriving this)
/// \brief Current energy of this bias (colvar_forces should be obtained by deriving this)
cvm::real bias_energy;
/// Whether to write the current bias energy from this bias to the trajectory file
bool b_output_energy;
/// \brief Whether this bias has already accumulated information
/// (when relevant)
bool has_data;
@ -94,6 +103,12 @@ public:
/// Write the bias configuration to a restart file
virtual std::ostream & write_restart (std::ostream &os);
/// Write a label to the trajectory file (comment line)
virtual std::ostream & write_traj_label (std::ostream &os);
/// Output quantities such as the bias energy to the trajectory file
virtual std::ostream & write_traj (std::ostream &os);
/// \brief Constructor
colvarbias_harmonic (std::string const &conf, char const *key);
@ -109,27 +124,9 @@ protected:
/// \brief Restraint centers without wrapping or constraints applied
std::vector<colvarvalue> colvar_centers_raw;
/// \brief Restraint force constant
cvm::real force_k;
/// \brief Moving target?
bool b_chg_centers;
/// \brief Changing force constant?
bool b_chg_force_k;
/// \brief Restraint force constant (target value)
cvm::real target_force_k;
/// \brief Equilibration steps for restraint FE calculation through TI
cvm::real target_equil_steps;
/// \brief Restraint force constant (starting value)
cvm::real starting_force_k;
/// \brief Lambda-schedule for custom varying force constant
std::vector<cvm::real> lambda_schedule;
/// \brief New restraint centers
std::vector<colvarvalue> target_centers;
@ -137,12 +134,41 @@ protected:
/// (or stage) towards the new values (calculated from target_nsteps)
std::vector<colvarvalue> centers_incr;
/// Whether to write the current restraint centers to the trajectory file
bool b_output_centers;
/// Whether to write the current accumulated work to the trajectory file
bool b_output_acc_work;
/// \brief Accumulated work
cvm::real acc_work;
/// \brief Restraint force constant
cvm::real force_k;
/// \brief Changing force constant?
bool b_chg_force_k;
/// \brief Restraint force constant (target value)
cvm::real target_force_k;
/// \brief Restraint force constant (starting value)
cvm::real starting_force_k;
/// \brief Lambda-schedule for custom varying force constant
std::vector<cvm::real> lambda_schedule;
/// \brief Exponent for varying the force constant
cvm::real force_k_exp;
/// \brief Intermediate quantity to compute the restraint free energy
/// (in TI, would be the accumulating FE derivative)
cvm::real restraint_FE;
/// \brief Number of steps required to reach the target force constant
/// or restraint centers
size_t target_nsteps;
/// \brief Equilibration steps for restraint FE calculation through TI
cvm::real target_equil_steps;
/// \brief Number of stages over which to perform the change
/// If zero, perform a continuous change
@ -150,10 +176,10 @@ protected:
/// \brief Number of current stage of the perturbation
int stage;
/// \brief Intermediate quantity to compute the restraint free energy
/// (in TI, would be the accumulating FE derivative)
cvm::real restraint_FE;
/// \brief Number of steps required to reach the target force constant
/// or restraint centers
size_t target_nsteps;
};

73
lib/colvars/colvarmodule.cpp
View File

@ -44,6 +44,11 @@ colvarmodule::colvarmodule (char const *config_filename,
config_s.close();
}
std::string index_file_name;
if (parse->get_keyval (conf, "indexFile", index_file_name)) {
read_index_file (index_file_name.c_str());
}
parse->get_keyval (conf, "analysis", b_analysis, false);
parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07,
@ -478,6 +483,11 @@ void colvarmodule::calc() {
cvi++) {
(*cvi)->write_traj_label (cv_traj_os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_traj_label (cv_traj_os);
}
cv_traj_os << "\n";
if (cvm::debug())
cv_traj_os.flush();
@ -494,6 +504,11 @@ void colvarmodule::calc() {
cvi++) {
(*cvi)->write_traj (cv_traj_os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_traj (cv_traj_os);
}
cv_traj_os << "\n";
if (cvm::debug())
cv_traj_os.flush();
@ -564,7 +579,6 @@ colvarmodule::~colvarmodule()
}
delete parse;
proxy = NULL;
}
@ -721,6 +735,55 @@ void cvm::exit (std::string const &message)
}
void cvm::read_index_file (char const *filename)
{
std::ifstream is (filename);
if (!is.good())
fatal_error ("Error: in opening index file \""+
std::string (filename)+"\".\n");
// std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
// std::list<std::vector<int> >::iterator lists_i = cvm::index_groups.begin();
while (is.good()) {
char open, close;
std::string group_name;
if ( (is >> open) && (open == '[') &&
(is >> group_name) &&
(is >> close) && (close == ']') ) {
cvm::index_group_names.push_back (group_name);
cvm::index_groups.push_back (std::vector<int> ());
} else {
cvm::fatal_error ("Error: in parsing index file \""+
std::string (filename)+"\".\n");
}
int atom_number = 1;
size_t pos = is.tellg();
while ( (is >> atom_number) && (atom_number > 0) ) {
(cvm::index_groups.back()).push_back (atom_number);
pos = is.tellg();
}
is.clear();
is.seekg (pos, std::ios::beg);
std::string delim;
if ( (is >> delim) && (delim == "[") ) {
// new group
is.clear();
is.seekg (pos, std::ios::beg);
} else {
break;
}
}
cvm::log ("The following index groups were read from the index file \""+
std::string (filename)+"\":\n");
std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
std::list<std::vector<int> >::iterator lists_i = cvm::index_groups.begin();
for ( ; names_i != cvm::index_group_names.end() ; names_i++, lists_i++) {
cvm::log (" "+(*names_i)+" ("+cvm::to_str (lists_i->size())+" atoms).\n");
}
}
// static pointers
std::vector<colvar *> colvarmodule::colvars;
@ -741,6 +804,8 @@ size_t colvarmodule::cv_traj_freq = 0;
size_t colvarmodule::depth = 0;
bool colvarmodule::b_analysis = false;
cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04;
std::list<std::string> colvarmodule::index_group_names;
std::list<std::vector<int> > colvarmodule::index_groups;
// file name prefixes
@ -932,15 +997,13 @@ cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
// Calculate the optimal rotation between two groups, and implement it
// as a quaternion. The method is the one documented in: Coutsias EA,
// Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput
// Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254
void colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2,
matrix2d<cvm::real, 4, 4> &S)

12
lib/colvars/colvarmodule.h
View File

@ -2,7 +2,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2013-04-23"
#define COLVARS_VERSION "2013-05-14"
#endif
#ifndef COLVARS_DEBUG
@ -290,6 +290,16 @@ public:
atom_pos const &ref_pos);
/// \brief Names of groups from a Gromacs .ndx file to be read at startup
static std::list<std::string> index_group_names;
/// \brief Groups from a Gromacs .ndx file read at startup
static std::list<std::vector<int> > index_groups;
/// \brief Read a Gromacs .ndx file
static void read_index_file (char const *filename);
/// \brief Create atoms from a file \param filename name of the file
/// (usually a PDB) \param atoms array of the atoms to be allocated
/// \param pdb_field (optiona) if "filename" is a PDB file, use this

2
lib/colvars/colvarvalue.cpp
View File

@ -242,7 +242,7 @@ std::istream & operator >> (std::istream &is, colvarvalue &x)
}
size_t colvarvalue::output_width (size_t const &real_width)
size_t colvarvalue::output_width (size_t const &real_width) const
{
switch (this->value_type) {
case colvarvalue::type_scalar:

4
lib/colvars/colvarvalue.h
View File

@ -265,10 +265,10 @@ public:
/// with a different type to this object
void error_rside (Type const &vt) const;
///<EFBFBD>Give the number of characters required to output this
/// Give the number of characters required to output this
/// colvarvalue, given the current type assigned and the number of
/// characters for a real number
size_t output_width (size_t const &real_width);
size_t output_width (size_t const &real_width) const;
// optimized routines for operations with an array; xv and inner as