diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp new file mode 100644 index 0000000000..f5e2ae8e76 --- /dev/null +++ b/lib/colvars/colvarbias_restraint.cpp @@ -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 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; +} diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h new file mode 100644 index 0000000000..1e55c6e079 --- /dev/null +++ b/lib/colvars/colvarbias_restraint.h @@ -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 colvar_centers; + + /// \brief Restraint centers without wrapping or constraints applied + std::vector colvar_centers_raw; + + /// \brief Moving target? + bool b_chg_centers; + + /// \brief New restraint centers + std::vector target_centers; + + /// \brief Amplitude of the restraint centers' increment at each step + /// (or stage) towards the new values (calculated from target_nsteps) + std::vector 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 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 diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp new file mode 100644 index 0000000000..5b6c333f7a --- /dev/null +++ b/lib/colvars/colvarscript.cpp @@ -0,0 +1,304 @@ +// -*- c++ -*- + +#include +#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])+" [args...]\n\ +\n\ +Managing the colvars module:\n\ + configfile -- read configuration from a file\n\ + config -- 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 -- 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 -- set frame number\n"; + } + + result += "\n\ +Accessing collective variables:\n\ + colvar value -- return the current value of the colvar \n\ + colvar update -- recalculate the colvar \n\ + colvar delete -- delete the colvar \n\ + colvar addforce -- apply given force on \n\ +\n\ +Accessing biases:\n\ + bias energy -- return the current energy of the bias \n\ + bias update -- recalculate the bias \n\ + bias delete -- delete the bias \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::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::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; +} diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h new file mode 100644 index 0000000000..4a9f37a2f7 --- /dev/null +++ b/lib/colvars/colvarscript.h @@ -0,0 +1,47 @@ +/// -*- c++ -*- + +#ifndef COLVARSCRIPT_H +#define COLVARSCRIPT_H + +#include +#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 diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp new file mode 100644 index 0000000000..53dc1c1c28 --- /dev/null +++ b/lib/colvars/colvartypes.cpp @@ -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 const &pos1, + std::vector const &pos2, + matrix2d &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 &S, + cvm::real S_eigval[4], + matrix2d &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 const &pos1, + std::vector const &pos2) +{ + matrix2d S; + matrix2d S_backup; + cvm::real S_eigval[4]; + matrix2d 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 &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 &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 &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 &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 S_new; + cvm::real S_new_eigval[4]; + matrix2d 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 b (n, 0.0); + std::vector z (n, 0.0); + + for (ip=0;ip 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= p) p=d[k=j]; + if (k != i) { + d[k]=d[i]; + d[i]=p; + for (j=0;j