jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12613 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2014-10-07 20:30:53 +00:00
parent 044084847b
commit 211c12c39e
5 changed files with 1649 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

530
lib/colvars/colvarbias_restraint.cpp Normal file
View File

@ -0,0 +1,530 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias_restraint.h"
colvarbias_restraint::colvarbias_restraint (std::string const &conf,
char const *key)
: colvarbias (conf, key),
target_nstages (0),
target_nsteps (0)
{
get_keyval (conf, "forceConstant", force_k, 1.0);
// get the initial restraint centers
colvar_centers.resize (colvars.size());
colvar_centers_raw.resize (colvars.size());
for (size_t i = 0; i < colvars.size(); i++) {
colvar_centers[i].type (colvars[i]->type());
colvar_centers_raw[i].type (colvars[i]->type());
}
if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
for (size_t i = 0; i < colvars.size(); i++) {
colvar_centers[i].apply_constraints();
colvar_centers_raw[i] = colvar_centers[i];
}
} else {
colvar_centers.clear();
cvm::error ("Error: must define the initial centers of the restraints.\n");
}
if (colvar_centers.size() != colvars.size())
cvm::error ("Error: number of centers does not match "
"that of collective variables.\n");
if (get_keyval (conf, "targetCenters", target_centers, colvar_centers)) {
b_chg_centers = true;
for (size_t i = 0; i < target_centers.size(); i++) {
target_centers[i].apply_constraints();
}
} else {
b_chg_centers = false;
target_centers.clear();
}
if (get_keyval (conf, "targetForceConstant", target_force_k, 0.0)) {
if (b_chg_centers)
cvm::error ("Error: cannot specify both targetCenters and targetForceConstant.\n");
starting_force_k = force_k;
b_chg_force_k = true;
get_keyval (conf, "targetEquilSteps", target_equil_steps, 0);
get_keyval (conf, "lambdaSchedule", lambda_schedule, lambda_schedule);
if (lambda_schedule.size()) {
// There is one more lambda-point than stages
target_nstages = lambda_schedule.size() - 1;
}
} else {
b_chg_force_k = false;
}
if (b_chg_centers || b_chg_force_k) {
get_keyval (conf, "targetNumSteps", target_nsteps, 0);
if (!target_nsteps)
cvm::error ("Error: targetNumSteps must be non-zero.\n");
if (get_keyval (conf, "targetNumStages", target_nstages, target_nstages) &&
lambda_schedule.size()) {
cvm::error ("Error: targetNumStages and lambdaSchedule are incompatible.\n");
}
if (target_nstages) {
// This means that either numStages of lambdaSchedule has been provided
stage = 0;
restraint_FE = 0.0;
}
if (get_keyval (conf, "targetForceExponent", force_k_exp, 1.0)) {
if (! b_chg_force_k)
cvm::log ("Warning: not changing force constant: targetForceExponent will be ignored\n");
if (force_k_exp < 1.0)
cvm::log ("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n");
}
}
get_keyval (conf, "outputCenters", b_output_centers, false);
if (b_chg_centers) {
get_keyval (conf, "outputAccumulatedWork", b_output_acc_work, false);
} else {
b_output_acc_work = false;
}
acc_work = 0.0;
if (cvm::debug())
cvm::log ("Done initializing a new restraint bias.\n");
}
colvarbias_restraint::~colvarbias_restraint ()
{
if (cvm::n_rest_biases > 0)
cvm::n_rest_biases -= 1;
}
void colvarbias_restraint::change_configuration (std::string const &conf)
{
get_keyval (conf, "forceConstant", force_k, force_k);
if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
for (size_t i = 0; i < colvars.size(); i++) {
colvar_centers[i].apply_constraints();
colvar_centers_raw[i] = colvar_centers[i];
}
}
}
cvm::real colvarbias_restraint::energy_difference (std::string const &conf)
{
std::vector<colvarvalue> alt_colvar_centers;
cvm::real alt_force_k;
cvm::real alt_bias_energy = 0.0;
get_keyval (conf, "forceConstant", alt_force_k, force_k);
alt_colvar_centers.resize (colvars.size());
size_t i;
for (i = 0; i < colvars.size(); i++) {
alt_colvar_centers[i].type (colvars[i]->type());
}
if (get_keyval (conf, "centers", alt_colvar_centers, colvar_centers)) {
for (i = 0; i < colvars.size(); i++) {
colvar_centers[i].apply_constraints();
}
}
for (i = 0; i < colvars.size(); i++) {
alt_bias_energy += restraint_potential(restraint_convert_k(alt_force_k, colvars[i]->width),
colvars[i],
alt_colvar_centers[i]);
}
return alt_bias_energy - bias_energy;
}
cvm::real colvarbias_restraint::update()
{
bias_energy = 0.0;
if (cvm::debug())
cvm::log ("Updating the restraint bias \""+this->name+"\".\n");
// Setup first stage of staged variable force constant calculation
if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) {
cvm::real lambda;
if (lambda_schedule.size()) {
lambda = lambda_schedule[0];
} else {
lambda = 0.0;
}
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow (lambda, force_k_exp);
cvm::log ("Restraint " + this->name + ", stage " +
cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda));
cvm::log ("Setting force constant to " + cvm::to_str (force_k));
}
if (b_chg_centers) {
if (!centers_incr.size()) {
// if this is the first calculation, calculate the advancement
// at each simulation step (or stage, if applicable)
// (take current stage into account: it can be non-zero
// if we are restarting a staged calculation)
centers_incr.resize (colvars.size());
for (size_t i = 0; i < colvars.size(); i++) {
centers_incr[i].type (colvars[i]->type());
centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) /
cvm::real ( target_nstages ? (target_nstages - stage) :
(target_nsteps - cvm::step_absolute()));
}
if (cvm::debug()) {
cvm::log ("Center increment for the restraint bias \""+
this->name+"\": "+cvm::to_str (centers_incr)+" at stage "+cvm::to_str (stage)+ ".\n");
}
}
if (target_nstages) {
if ((cvm::step_relative() > 0)
&& (cvm::step_absolute() % target_nsteps) == 0
&& stage < target_nstages) {
for (size_t i = 0; i < colvars.size(); i++) {
colvar_centers_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[i];
colvars[i]->wrap(colvar_centers[i]);
colvar_centers[i].apply_constraints();
}
stage++;
cvm::log ("Moving restraint \"" + this->name +
"\" stage " + cvm::to_str(stage) +
" : setting centers to " + cvm::to_str (colvar_centers) +
" at step " + cvm::to_str (cvm::step_absolute()));
}
} else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) {
// move the restraint centers in the direction of the targets
// (slow growth)
for (size_t i = 0; i < colvars.size(); i++) {
colvar_centers_raw[i] += centers_incr[i];
colvar_centers[i] = colvar_centers_raw[i];
colvars[i]->wrap(colvar_centers[i]);
colvar_centers[i].apply_constraints();
}
}
if (cvm::debug())
cvm::log ("Current centers for the restraint bias \""+
this->name+"\": "+cvm::to_str (colvar_centers)+".\n");
}
if (b_chg_force_k) {
// Coupling parameter, between 0 and 1
cvm::real lambda;
if (target_nstages) {
// TI calculation: estimate free energy derivative
// need current lambda
if (lambda_schedule.size()) {
lambda = lambda_schedule[stage];
} else {
lambda = cvm::real(stage) / cvm::real(target_nstages);
}
if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) {
// Start averaging after equilibration period, if requested
// Square distance normalized by square colvar width
cvm::real dist_sq = 0.0;
for (size_t i = 0; i < colvars.size(); i++) {
dist_sq += colvars[i]->dist2 (colvars[i]->value(), colvar_centers[i])
/ (colvars[i]->width * colvars[i]->width);
}
restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0)
* (target_force_k - starting_force_k) * dist_sq;
}
// Finish current stage...
if (cvm::step_absolute() % target_nsteps == 0 &&
cvm::step_absolute() > 0) {
cvm::log ("Lambda= " + cvm::to_str (lambda) + " dA/dLambda= "
+ cvm::to_str (restraint_FE / cvm::real(target_nsteps - target_equil_steps)));
// ...and move on to the next one
if (stage < target_nstages) {
restraint_FE = 0.0;
stage++;
if (lambda_schedule.size()) {
lambda = lambda_schedule[stage];
} else {
lambda = cvm::real(stage) / cvm::real(target_nstages);
}
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow (lambda, force_k_exp);
cvm::log ("Restraint " + this->name + ", stage " +
cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda));
cvm::log ("Setting force constant to " + cvm::to_str (force_k));
}
}
} else if (cvm::step_absolute() <= target_nsteps) {
// update force constant (slow growth)
lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps);
force_k = starting_force_k + (target_force_k - starting_force_k)
* std::pow (lambda, force_k_exp);
}
}
if (cvm::debug())
cvm::log ("Done updating the restraint bias \""+this->name+"\".\n");
// Force and energy calculation
for (size_t i = 0; i < colvars.size(); i++) {
colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(force_k, colvars[i]->width),
colvars[i],
colvar_centers[i]);
bias_energy += restraint_potential(restraint_convert_k(force_k, colvars[i]->width),
colvars[i],
colvar_centers[i]);
if (cvm::debug()) {
cvm::log ("dist_grad["+cvm::to_str (i)+
"] = "+cvm::to_str (colvars[i]->dist2_lgrad (colvars[i]->value(),
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 restraint bias \""+
this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
return bias_energy;
}
std::istream & colvarbias_restraint::read_restart (std::istream &is)
{
size_t const start_pos = is.tellg();
cvm::log ("Restarting restraint bias \""+
this->name+"\".\n");
std::string key, brace, conf;
if ( !(is >> key) || !(key == "restraint" || key == "harmonic") ||
!(is >> brace) || !(brace == "{") ||
!(is >> colvarparse::read_block ("configuration", conf)) ) {
cvm::log ("Error: in reading restart configuration for restraint bias \""+
this->name+"\" at position "+
cvm::to_str (is.tellg())+" in stream.\n");
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
// int id = -1;
std::string name = "";
// if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) &&
// (id != this->id) ) ||
if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) &&
(name != this->name) )
cvm::error ("Error: in the restart file, the "
"\"restraint\" block has a wrong name\n");
// if ( (id == -1) && (name == "") ) {
if (name.size() == 0) {
cvm::error ("Error: \"restraint\" block in the restart file "
"has no identifiers.\n");
}
if (b_chg_centers) {
// cvm::log ("Reading the updated restraint centers from the restart.\n");
if (!get_keyval (conf, "centers", colvar_centers))
cvm::error ("Error: restraint centers are missing from the restart.\n");
if (!get_keyval (conf, "centers_raw", colvar_centers_raw))
cvm::error ("Error: \"raw\" restraint centers are missing from the restart.\n");
}
if (b_chg_force_k) {
// cvm::log ("Reading the updated force constant from the restart.\n");
if (!get_keyval (conf, "forceConstant", force_k))
cvm::error ("Error: force constant is missing from the restart.\n");
}
if (target_nstages) {
// cvm::log ("Reading current stage from the restart.\n");
if (!get_keyval (conf, "stage", stage))
cvm::error ("Error: current stage is missing from the restart.\n");
}
if (b_output_acc_work) {
if (!get_keyval (conf, "accumulatedWork", acc_work))
cvm::error ("Error: accumulatedWork is missing from the restart.\n");
}
is >> brace;
if (brace != "}") {
cvm::error ("Error: corrupt restart information for restraint bias \""+
this->name+"\": no matching brace at position "+
cvm::to_str (is.tellg())+" in the restart file.\n");
is.setstate (std::ios::failbit);
}
return is;
}
std::ostream & colvarbias_restraint::write_restart (std::ostream &os)
{
os << "restraint {\n"
<< " configuration {\n"
// << " id " << this->id << "\n"
<< " name " << this->name << "\n";
if (b_chg_centers) {
size_t i;
os << " centers ";
for (i = 0; i < colvars.size(); i++) {
os << " " << colvar_centers[i];
}
os << "\n";
os << " centers_raw ";
for (i = 0; i < colvars.size(); i++) {
os << " " << colvar_centers_raw[i];
}
os << "\n";
}
if (b_chg_force_k) {
os << " forceConstant "
<< std::setprecision (cvm::en_prec)
<< std::setw (cvm::en_width) << force_k << "\n";
}
if (target_nstages) {
os << " stage " << std::setw (cvm::it_width)
<< stage << "\n";
}
if (b_output_acc_work) {
os << " accumulatedWork " << acc_work << "\n";
}
os << " }\n"
<< "}\n\n";
return os;
}
std::ostream & colvarbias_restraint::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_restraint::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;
}
colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(std::string const &conf, char const *key) :
colvarbias_restraint(conf, key) {
for (size_t i = 0; i < colvars.size(); i++) {
if (colvars[i]->width != 1.0)
cvm::log ("The force constant for colvar \""+colvars[i]->name+
"\" will be rescaled to "+
cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+
" according to the specified width.\n");
}
}
cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const
{
return 0.5 * k * x->dist2(x->value(), xcenter);
}
colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const
{
return 0.5 * k * x->dist2_lgrad(x->value(), xcenter);
}
cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
{
return k / (dist_measure * dist_measure);
}
colvarbias_restraint_linear::colvarbias_restraint_linear(std::string const &conf, char const *key) :
colvarbias_restraint(conf, key) {
for (size_t i = 0; i < colvars.size(); i++) {
if (colvars[i]->width != 1.0)
cvm::log ("The force constant for colvar \""+colvars[i]->name+
"\" will be rescaled to "+
cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+
" according to the specified width.\n");
}
}
cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const
{
return k * (x->value() - xcenter);
}
colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const
{
return k;
}
cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, cvm::real dist_measure) const
{
return k / dist_measure;
}

153
lib/colvars/colvarbias_restraint.h Normal file
View File

@ -0,0 +1,153 @@
/// -*- c++ -*-
#ifndef COLVARBIAS_RESTRAINT_H
#define COLVARBIAS_RESTRAINT_H
#include "colvarbias.h"
/// \brief Bias restraint, optionally moving towards a target
/// (implementation of \link colvarbias \endlink)
class colvarbias_restraint : public colvarbias {
public:
/// Retrieve colvar values and calculate their biasing forces
virtual cvm::real update();
/// Load new configuration - force constant and/or centers only
virtual void change_configuration(std::string const &conf);
/// Calculate change in energy from using alternate configuration
virtual cvm::real energy_difference(std::string const &conf);
/// Read the bias configuration from a restart file
virtual std::istream & read_restart (std::istream &is);
/// 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_restraint (std::string const &conf, char const *key);
/// Destructor
virtual ~colvarbias_restraint();
protected:
/// \brief Potential function
virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0;
/// \brief Force function
virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0;
///\brief Unit scaling
virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0;
/// \brief Restraint centers
std::vector<colvarvalue> colvar_centers;
/// \brief Restraint centers without wrapping or constraints applied
std::vector<colvarvalue> colvar_centers_raw;
/// \brief Moving target?
bool b_chg_centers;
/// \brief New restraint centers
std::vector<colvarvalue> target_centers;
/// \brief Amplitude of the restraint centers' increment at each step
/// (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 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
int target_nstages;
/// \brief Number of current stage of the perturbation
int stage;
/// \brief Number of steps required to reach the target force constant
/// or restraint centers
size_t target_nsteps;
};
/// \brief Harmonic bias restraint
/// (implementation of \link colvarbias_restraint \endlink)
class colvarbias_restraint_harmonic : public colvarbias_restraint {
public:
colvarbias_restraint_harmonic(std::string const &conf, char const *key);
protected: /// \brief Potential function
virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const;
/// \brief Force function
virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const;
///\brief Unit scaling
virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
};
/// \brief Linear bias restraint
/// (implementation of \link colvarbias_restraint \endlink)
class colvarbias_restraint_linear : public colvarbias_restraint {
public:
colvarbias_restraint_linear(std::string const &conf, char const *key);
protected: /// \brief Potential function
virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const;
/// \brief Force function
virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const;
///\brief Unit scaling
virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const;
};
#endif

304
lib/colvars/colvarscript.cpp Normal file
View File

@ -0,0 +1,304 @@
// -*- c++ -*-
#include <cstdlib>
#include "colvarscript.h"
colvarscript::colvarscript (colvarproxy *p)
: proxy (p),
colvars (p->colvars),
proxy_error (0)
{
}
/// Run method based on given arguments
int colvarscript::run (int argc, char const *argv[]) {
result = "";
if (cvm::debug()) {
cvm::log ("Called script run with " + cvm::to_str(argc) + " args");
for (int i = 0; i < argc; i++) { cvm::log (argv[i]); }
}
if (argc < 2) {
result = "usage: "+std::string (argv[0])+" <subcommand> [args...]\n\
\n\
Managing the colvars module:\n\
configfile <file name> -- read configuration from a file\n\
config <string> -- read configuration from the given string\n\
reset -- delete all internal configuration\n\
delete -- delete this colvars module instance\n\
\n\
Input and output:\n\
list -- return a list of all variables\n\
list biases -- return a list of all biases\n\
load <file name> -- load a state file (requires configuration)\n\
update -- recalculate colvars and biases based\n\
printframe -- return a summary of the current frame\n\
printframelabels -- return labels to annotate printframe's output\n";
if (proxy->frame() != COLVARS_NOT_IMPLEMENTED) {
result += "\
frame -- return current frame number\n\
frame <new_frame> -- set frame number\n";
}
result += "\n\
Accessing collective variables:\n\
colvar <name> value -- return the current value of the colvar <name>\n\
colvar <name> update -- recalculate the colvar <name>\n\
colvar <name> delete -- delete the colvar <name>\n\
colvar <name> addforce <F> -- apply given force on <name>\n\
\n\
Accessing biases:\n\
bias <name> energy -- return the current energy of the bias <name>\n\
bias <name> update -- recalculate the bias <name>\n\
bias <name> delete -- delete the bias <name>\n\
\n\
";
return COLVARSCRIPT_OK;
}
std::string cmd = argv[1];
if (cmd == "colvar") {
return proc_colvar (argc-1, &(argv[1]));
}
if (cmd == "bias") {
return proc_bias (argc-1, &(argv[1]));
}
if (cmd == "reset") {
/// Delete every child object
colvars->reset();
return COLVARSCRIPT_OK;
}
if (cmd == "delete") {
colvars->reset();
// Note: the delete bit may be ignored by some backends
// it is mostly useful in VMD
colvars->set_error_bits(DELETE_COLVARS);
return COLVARSCRIPT_OK;
}
if (cmd == "update") {
colvars->calc();
return COLVARSCRIPT_OK;
}
if (cmd == "list") {
if (argc == 2) {
for (std::vector<colvar *>::iterator cvi = colvars->colvars.begin();
cvi != colvars->colvars.end();
++cvi) {
result += (cvi == colvars->colvars.begin() ? "" : " ") + (*cvi)->name;
}
return COLVARSCRIPT_OK;
} else if (argc == 3 && !strcmp(argv[2], "biases")) {
for (std::vector<colvarbias *>::iterator bi = colvars->biases.begin();
bi != colvars->biases.end();
++bi) {
result += (bi == colvars->biases.begin() ? "" : " ") + (*bi)->name;
}
return COLVARSCRIPT_OK;
} else {
result = "Wrong arguments to command \"list\"";
return COLVARSCRIPT_ERROR;
}
}
/// Parse config from file
if (cmd == "configfile") {
if (argc < 3) {
result = "Missing arguments";
return COLVARSCRIPT_ERROR;
}
if (colvars->config_file (argv[2]) == COLVARS_OK) {
return COLVARSCRIPT_OK;
} else {
return COLVARSCRIPT_ERROR;
}
}
/// Parse config from string
if (cmd == "config") {
if (argc < 3) {
result = "Missing arguments";
return COLVARSCRIPT_ERROR;
}
std::string conf = argv[2];
if (colvars->config_string (conf) == COLVARS_OK) {
return COLVARSCRIPT_OK;
} else {
return COLVARSCRIPT_ERROR;
}
}
/// Load an input state file
if (cmd == "load") {
if (argc < 3) {
result = "Missing arguments";
return COLVARSCRIPT_ERROR;
}
proxy->input_prefix_str = argv[2];
if (colvars->setup_input() == COLVARS_OK) {
return COLVARSCRIPT_OK;
} else {
return COLVARSCRIPT_ERROR;
}
}
/// TODO Write an output state file? (Useful for testing)
/// Print the values that would go on colvars.traj
if (cmd == "printframelabels") {
std::ostringstream os;
colvars->write_traj_label (os);
result = os.str();
return COLVARSCRIPT_OK;
}
if (cmd == "printframe") {
std::ostringstream os;
colvars->write_traj (os);
result = os.str();
return COLVARSCRIPT_OK;
}
if (cmd == "frame") {
if (argc == 2) {
int f = proxy->frame();
if (f >= 0) {
result = cvm::to_str (f);
return COLVARSCRIPT_OK;
} else {
result = "Frame number is not available";
return COLVARSCRIPT_ERROR;
}
} else if (argc == 3) {
// Failure of this function does not trigger an error, but
// returns the plain result to let scripts detect available frames
long int f = proxy->frame(strtol(argv[2], NULL, 10));
colvars->it = proxy->frame();
result = cvm::to_str (f);
return COLVARSCRIPT_OK;
} else {
result = "Wrong arguments to command \"frame\"";
return COLVARSCRIPT_ERROR;
}
}
result = "Syntax error";
return COLVARSCRIPT_ERROR;
}
int colvarscript::proc_colvar (int argc, char const *argv[]) {
std::string name = argv[1];
colvar *cv = cvm::colvar_by_name (name);
if (cv == NULL) {
result = "Colvar not found: " + name;
return COLVARSCRIPT_ERROR;
}
if (argc < 3) {
result = "Missing parameters";
return COLVARSCRIPT_ERROR;
}
std::string subcmd = argv[2];
if (subcmd == "value") {
result = cvm::to_str(cv->value(), 0, cvm::cv_prec);
return COLVARSCRIPT_OK;
}
if (subcmd == "update") {
// note: this is not sufficient for a newly created colvar
// as atom coordinates may not be properly loaded
// a full CVM update is required
// or otherwise updating atom positions
cv->update();
result = cvm::to_str(cv->value(), 0, cvm::cv_prec);
return COLVARSCRIPT_OK;
}
if (subcmd == "delete") {
if (cv->biases.size() > 0) {
result = "Cannot delete a colvar currently used by biases, delete those biases first";
return COLVARSCRIPT_ERROR;
}
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
colvars->write_traj_label (colvars->cv_traj_os);
return COLVARSCRIPT_OK;
}
if (subcmd == "addforce") {
if (argc < 4) {
result = "Missing parameter: force value";
return COLVARSCRIPT_ERROR;
}
std::string f_str = argv[3];
std::istringstream is (f_str);
is.width(cvm::cv_width);
is.precision(cvm::cv_prec);
colvarvalue force (cv->type());
force.is_derivative();
if (!(is >> force)) {
result = "Error parsing force value";
return COLVARSCRIPT_ERROR;
}
cv->add_bias_force(force);
result = cvm::to_str(force, cvm::cv_width, cvm::cv_prec);
return COLVARSCRIPT_OK;
}
result = "Syntax error";
return COLVARSCRIPT_ERROR;
}
int colvarscript::proc_bias (int argc, char const *argv[]) {
std::string name = argv[1];
colvarbias *b = cvm::bias_by_name (name);
if (b == NULL) {
result = "Bias not found: " + name;
return COLVARSCRIPT_ERROR;
}
if (argc < 3) {
result = "Missing parameters";
return COLVARSCRIPT_ERROR;
}
std::string subcmd = argv[2];
if (subcmd == "energy") {
result = cvm::to_str(b->get_energy());
return COLVARSCRIPT_OK;
}
if (subcmd == "update") {
b->update();
result = cvm::to_str(b->get_energy());
return COLVARSCRIPT_OK;
}
if (subcmd == "delete") {
// the bias destructor takes care of the cleanup at cvm level
delete b;
// TODO this could be done by the destructors
colvars->write_traj_label (colvars->cv_traj_os);
return COLVARSCRIPT_OK;
}
if (argc >= 4) {
// std::string param = argv[3];
result = "Syntax error";
return COLVARSCRIPT_ERROR;
}
result = "Syntax error";
return COLVARSCRIPT_ERROR;
}

47
lib/colvars/colvarscript.h Normal file
View File

@ -0,0 +1,47 @@
/// -*- c++ -*-
#ifndef COLVARSCRIPT_H
#define COLVARSCRIPT_H
#include <string>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarbias.h"
#include "colvarproxy.h"
#define COLVARSCRIPT_ERROR -1
#define COLVARSCRIPT_OK 0
class colvarscript {
private:
colvarproxy *proxy;
colvarmodule *colvars;
inline colvarscript() {} // no-argument construction forbidden
public:
friend class colvarproxy;
colvarscript(colvarproxy * p);
inline ~colvarscript() {}
/// If an error is caught by the proxy through fatal_error(), this is set to COLVARSCRIPT_ERROR
int proxy_error;
/// If an error is returned by one of the methods, it should set this to the error message
std::string result;
/// Run script command with given positional arguments
int run (int argc, char const *argv[]);
/// Run subcommands on colvar
int proc_colvar (int argc, char const *argv[]);
/// Run subcommands on bias
int proc_bias (int argc, char const *argv[]);
};
#endif

615
lib/colvars/colvartypes.cpp Normal file
View File

@ -0,0 +1,615 @@
/// -*- c++ -*-
#include "colvarmodule.h"
#include "colvartypes.h"
#include "colvarparse.h"
std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os << "( ";
os.width (w); os.precision (p);
os << v.x << " , ";
os.width (w); os.precision (p);
os << v.y << " , ";
os.width (w); os.precision (p);
os << v.z << " )";
return os;
}
std::istream & operator >> (std::istream &is, colvarmodule::rvector &v)
{
size_t const start_pos = is.tellg();
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> v.x) || !(is >> sep) || !(sep == ',') ||
!(is >> v.y) || !(is >> sep) || !(sep == ',') ||
!(is >> v.z) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
return is;
}
std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os << "( ";
os.width (w); os.precision (p);
os << q.q0 << " , ";
os.width (w); os.precision (p);
os << q.q1 << " , ";
os.width (w); os.precision (p);
os << q.q2 << " , ";
os.width (w); os.precision (p);
os << q.q3 << " )";
return os;
}
std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
{
size_t const start_pos = is.tellg();
std::string euler ("");
if ( (is >> euler) && (colvarparse::to_lower_cppstr (euler) ==
std::string ("euler")) ) {
// parse the Euler angles
char sep;
cvm::real phi, theta, psi;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> phi) || !(is >> sep) || !(sep == ',') ||
!(is >> theta) || !(is >> sep) || !(sep == ',') ||
!(is >> psi) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
q = colvarmodule::quaternion (phi, theta, psi);
} else {
// parse the quaternion components
is.seekg (start_pos, std::ios::beg);
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> q.q0) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q1) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q2) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q3) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
}
return is;
}
cvm::quaternion
cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
cvm::rvector const &vec) const
{
cvm::quaternion result (0.0, 0.0, 0.0, 0.0);
result.q0 = 2.0 * pos.x * q0 * vec.x
+2.0 * pos.y * q0 * vec.y
+2.0 * pos.z * q0 * vec.z
-2.0 * pos.y * q3 * vec.x
+2.0 * pos.z * q2 * vec.x
+2.0 * pos.x * q3 * vec.y
-2.0 * pos.z * q1 * vec.y
-2.0 * pos.x * q2 * vec.z
+2.0 * pos.y * q1 * vec.z;
result.q1 = +2.0 * pos.x * q1 * vec.x
-2.0 * pos.y * q1 * vec.y
-2.0 * pos.z * q1 * vec.z
+2.0 * pos.y * q2 * vec.x
+2.0 * pos.z * q3 * vec.x
+2.0 * pos.x * q2 * vec.y
-2.0 * pos.z * q0 * vec.y
+2.0 * pos.x * q3 * vec.z
+2.0 * pos.y * q0 * vec.z;
result.q2 = -2.0 * pos.x * q2 * vec.x
+2.0 * pos.y * q2 * vec.y
-2.0 * pos.z * q2 * vec.z
+2.0 * pos.y * q1 * vec.x
+2.0 * pos.z * q0 * vec.x
+2.0 * pos.x * q1 * vec.y
+2.0 * pos.z * q3 * vec.y
-2.0 * pos.x * q0 * vec.z
+2.0 * pos.y * q3 * vec.z;
result.q3 = -2.0 * pos.x * q3 * vec.x
-2.0 * pos.y * q3 * vec.y
+2.0 * pos.z * q3 * vec.z
-2.0 * pos.y * q0 * vec.x
+2.0 * pos.z * q1 * vec.x
+2.0 * pos.x * q0 * vec.y
+2.0 * pos.z * q2 * vec.y
+2.0 * pos.x * q1 * vec.z
+2.0 * pos.y * q2 * vec.z;
return result;
}
// Calculate the optimal rotation between two groups, and implement it
// as a quaternion. Uses the method 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)
{
cvm::rmatrix C;
// build the correlation matrix
C.reset();
for (size_t i = 0; i < pos1.size(); i++) {
C.xx() += pos1[i].x * pos2[i].x;
C.xy() += pos1[i].x * pos2[i].y;
C.xz() += pos1[i].x * pos2[i].z;
C.yx() += pos1[i].y * pos2[i].x;
C.yy() += pos1[i].y * pos2[i].y;
C.yz() += pos1[i].y * pos2[i].z;
C.zx() += pos1[i].z * pos2[i].x;
C.zy() += pos1[i].z * pos2[i].y;
C.zz() += pos1[i].z * pos2[i].z;
}
// build the "overlap" matrix, whose eigenvectors are stationary
// points of the RMSD in the space of rotations
S[0][0] = C.xx() + C.yy() + C.zz();
S[1][0] = C.yz() - C.zy();
S[0][1] = S[1][0];
S[2][0] = - C.xz() + C.zx() ;
S[0][2] = S[2][0];
S[3][0] = C.xy() - C.yx();
S[0][3] = S[3][0];
S[1][1] = C.xx() - C.yy() - C.zz();
S[2][1] = C.xy() + C.yx();
S[1][2] = S[2][1];
S[3][1] = C.xz() + C.zx();
S[1][3] = S[3][1];
S[2][2] = - C.xx() + C.yy() - C.zz();
S[3][2] = C.yz() + C.zy();
S[2][3] = S[3][2];
S[3][3] = - C.xx() - C.yy() + C.zz();
// if (cvm::debug()) {
// for (size_t i = 0; i < 4; i++) {
// std::string line ("");
// for (size_t j = 0; j < 4; j++) {
// line += std::string (" S["+cvm::to_str (i)+
// "]["+cvm::to_str (j)+"] ="+cvm::to_str (S[i][j]));
// }
// cvm::log (line+"\n");
// }
// }
}
void colvarmodule::rotation::diagonalize_matrix (matrix2d<cvm::real, 4, 4> &S,
cvm::real S_eigval[4],
matrix2d<cvm::real, 4, 4> &S_eigvec)
{
// diagonalize
int jac_nrot = 0;
jacobi (S, S_eigval, S_eigvec, &jac_nrot);
eigsrt (S_eigval, S_eigvec);
// jacobi saves eigenvectors by columns
transpose (S_eigvec);
// normalize eigenvectors
for (size_t ie = 0; ie < 4; ie++) {
cvm::real norm2 = 0.0;
size_t i;
for (i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2));
cvm::real const norm = std::sqrt (norm2);
for (i = 0; i < 4; i++) S_eigvec[ie][i] /= norm;
}
}
// Calculate the rotation, plus its derivatives
void colvarmodule::rotation::calc_optimal_rotation
(std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2)
{
matrix2d<cvm::real, 4, 4> S;
matrix2d<cvm::real, 4, 4> S_backup;
cvm::real S_eigval[4];
matrix2d<cvm::real, 4, 4> S_eigvec;
// if (cvm::debug()) {
// cvm::atom_pos cog1 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos1.size(); i++) {
// cog1 += pos1[i];
// }
// cog1 /= cvm::real (pos1.size());
// cvm::atom_pos cog2 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos2.size(); i++) {
// cog2 += pos2[i];
// }
// cog2 /= cvm::real (pos1.size());
// cvm::log ("calc_optimal_rotation: centers of geometry are: "+
// cvm::to_str (cog1, cvm::cv_width, cvm::cv_prec)+
// " and "+cvm::to_str (cog2, cvm::cv_width, cvm::cv_prec)+".\n");
// }
build_matrix (pos1, pos2, S);
S_backup = S;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("S = "+cvm::to_str (cvm::to_str (S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
}
}
diagonalize_matrix (S, S_eigval, S_eigvec);
// eigenvalues and eigenvectors
cvm::real const &L0 = S_eigval[0];
cvm::real const &L1 = S_eigval[1];
cvm::real const &L2 = S_eigval[2];
cvm::real const &L3 = S_eigval[3];
cvm::real const *Q0 = S_eigvec[0];
cvm::real const *Q1 = S_eigvec[1];
cvm::real const *Q2 = S_eigvec[2];
cvm::real const *Q3 = S_eigvec[3];
lambda = L0;
q = cvm::quaternion (Q0);
if (q_old.norm2() > 0.0) {
q.match (q_old);
if (q_old.inner (q) < (1.0 - crossing_threshold)) {
cvm::log ("Warning: one molecular orientation has changed by more than "+
cvm::to_str (crossing_threshold)+": discontinuous rotation ?\n");
}
}
q_old = q;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str (cvm::quaternion (Q0), cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q0)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L1 = "+cvm::to_str (L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str (cvm::quaternion (Q1), cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q1)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L2 = "+cvm::to_str (L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str (cvm::quaternion (Q2), cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q2)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L3 = "+cvm::to_str (L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str (cvm::quaternion (Q3), cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q3)), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
// calculate derivatives of L0 and Q0 with respect to each atom in
// either group; note: if dS_1 is a null vector, nothing will be
// calculated
size_t ia;
for (ia = 0; ia < dS_1.size(); ia++) {
cvm::real const &a2x = pos2[ia].x;
cvm::real const &a2y = pos2[ia].y;
cvm::real const &a2z = pos2[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_1 = dS_1[ia];
// derivative of the S matrix
ds_1.reset();
ds_1[0][0] = cvm::rvector ( a2x, a2y, a2z);
ds_1[1][0] = cvm::rvector ( 0.0, a2z, -a2y);
ds_1[0][1] = ds_1[1][0];
ds_1[2][0] = cvm::rvector (-a2z, 0.0, a2x);
ds_1[0][2] = ds_1[2][0];
ds_1[3][0] = cvm::rvector ( a2y, -a2x, 0.0);
ds_1[0][3] = ds_1[3][0];
ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z);
ds_1[2][1] = cvm::rvector ( a2y, a2x, 0.0);
ds_1[1][2] = ds_1[2][1];
ds_1[3][1] = cvm::rvector ( a2z, 0.0, a2x);
ds_1[1][3] = ds_1[3][1];
ds_1[2][2] = cvm::rvector (-a2x, a2y, -a2z);
ds_1[3][2] = cvm::rvector ( 0.0, a2z, a2y);
ds_1[2][3] = ds_1[3][2];
ds_1[3][3] = cvm::rvector (-a2x, -a2y, a2z);
cvm::rvector &dl0_1 = dL0_1[ia];
vector1d<cvm::rvector, 4> &dq0_1 = dQ0_1[ia];
// matrix multiplications; derivatives of L_0 and Q_0 are
// calculated using Hellmann-Feynman theorem (i.e. exploiting the
// fact that the eigenvectors Q_i form an orthonormal basis)
dl0_1.reset();
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dl0_1 += Q0[i] * ds_1[i][j] * Q0[j];
}
}
dq0_1.reset();
for (size_t p = 0; p < 4; p++) {
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dq0_1[p] +=
(Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
(Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
(Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
}
}
}
}
// do the same for the second group
for (ia = 0; ia < dS_2.size(); ia++) {
cvm::real const &a1x = pos1[ia].x;
cvm::real const &a1y = pos1[ia].y;
cvm::real const &a1z = pos1[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_2 = dS_2[ia];
ds_2.reset();
ds_2[0][0] = cvm::rvector ( a1x, a1y, a1z);
ds_2[1][0] = cvm::rvector ( 0.0, -a1z, a1y);
ds_2[0][1] = ds_2[1][0];
ds_2[2][0] = cvm::rvector ( a1z, 0.0, -a1x);
ds_2[0][2] = ds_2[2][0];
ds_2[3][0] = cvm::rvector (-a1y, a1x, 0.0);
ds_2[0][3] = ds_2[3][0];
ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z);
ds_2[2][1] = cvm::rvector ( a1y, a1x, 0.0);
ds_2[1][2] = ds_2[2][1];
ds_2[3][1] = cvm::rvector ( a1z, 0.0, a1x);
ds_2[1][3] = ds_2[3][1];
ds_2[2][2] = cvm::rvector (-a1x, a1y, -a1z);
ds_2[3][2] = cvm::rvector ( 0.0, a1z, a1y);
ds_2[2][3] = ds_2[3][2];
ds_2[3][3] = cvm::rvector (-a1x, -a1y, a1z);
cvm::rvector &dl0_2 = dL0_2[ia];
vector1d<cvm::rvector, 4> &dq0_2 = dQ0_2[ia];
dl0_2.reset();
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dl0_2 += Q0[i] * ds_2[i][j] * Q0[j];
}
}
dq0_2.reset();
for (size_t p = 0; p < 4; p++) {
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dq0_2[p] +=
(Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
(Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
(Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p];
}
}
}
if (cvm::debug()) {
if (b_debug_gradients) {
matrix2d<cvm::real, 4, 4> S_new;
cvm::real S_new_eigval[4];
matrix2d<cvm::real, 4, 4> S_new_eigvec;
// make an infitesimal move along each cartesian coordinate of
// this atom, and solve again the eigenvector problem
for (size_t comp = 0; comp < 3; comp++) {
S_new = S_backup;
// diagonalize the new overlap matrix
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
S_new[i][j] +=
colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
}
}
// cvm::log ("S_new = "+cvm::to_str (cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix (S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::real const *Q0_new = S_new_eigvec[0];
cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
cvm::quaternion const q0 (Q0);
cvm::quaternion const DQ0 (dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::log ( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str (std::fabs (L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str ((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
}
}
}
// Numerical Recipes routine for diagonalization
#define ROTATE(a,i,j,k,l) g=a[i][j]; \
h=a[k][l]; \
a[i][j]=g-s*(h+g*tau); \
a[k][l]=h+s*(g-h*tau);
#define n 4
void jacobi(cvm::real **a, cvm::real d[], cvm::real **v, int *nrot)
{
int j,iq,ip,i;
cvm::real tresh,theta,tau,t,sm,s,h,g,c;
std::vector<cvm::real> b (n, 0.0);
std::vector<cvm::real> z (n, 0.0);
for (ip=0;ip<n;ip++) {
for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++) {
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=0;i<=50;i++) {
sm=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++)
sm += std::fabs(a[ip][iq]);
}
if (sm == 0.0) {
return;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++) {
g=100.0*std::fabs(a[ip][iq]);
if (i > 4 && (cvm::real)(std::fabs(d[ip])+g) == (cvm::real)std::fabs(d[ip])
&& (cvm::real)(std::fabs(d[iq])+g) == (cvm::real)std::fabs(d[iq]))
a[ip][iq]=0.0;
else if (std::fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h))
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
t=1.0/(std::fabs(theta)+std::sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/std::sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=0;j<=ip-1;j++) {
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<=iq-1;j++) {
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<n;j++) {
ROTATE(a,ip,j,iq,j)
}
for (j=0;j<n;j++) {
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=0;ip<n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
cvm::error ("Too many iterations in routine jacobi.\n");
}
void eigsrt(cvm::real d[], cvm::real **v)
{
int k,j,i;
cvm::real p;
for (i=0;i<n;i++) {
p=d[k=i];
for (j=i+1;j<n;j++)
if (d[j] >= p) p=d[k=j];
if (k != i) {
d[k]=d[i];
d[i]=p;
for (j=0;j<n;j++) {
p=v[j][i];
v[j][i]=v[j][k];
v[j][k]=p;
}
}
}
}
void transpose(cvm::real **v)
{
cvm::real p;
for (int i=0;i<n;i++) {
for (int j=i+1;j<n;j++) {
p=v[i][j];
v[i][j]=v[j][i];
v[j][i]=p;
}
}
}
#undef n
#undef ROTATE