forked from lijiext/lammps
1370 lines
39 KiB
C++
1370 lines
39 KiB
C++
#include "colvarmodule.h"
|
|
#include "colvarparse.h"
|
|
#include "colvarproxy.h"
|
|
#include "colvar.h"
|
|
#include "colvarbias.h"
|
|
#include "colvarbias_meta.h"
|
|
#include "colvarbias_abf.h"
|
|
|
|
|
|
colvarmodule::colvarmodule (char const *config_filename,
|
|
colvarproxy *proxy_in)
|
|
{
|
|
// pointer to the proxy object
|
|
if (proxy == NULL) {
|
|
proxy = proxy_in;
|
|
parse = new colvarparse();
|
|
} else {
|
|
cvm::fatal_error ("Error: trying to allocate twice the collective "
|
|
"variable module.\n");
|
|
}
|
|
|
|
cvm::log (cvm::line_marker);
|
|
cvm::log ("Initializing the collective variables module, version "+
|
|
cvm::to_str(COLVARS_VERSION)+".\n");
|
|
|
|
// "it_restart" will be set by the input state file, if any;
|
|
// "it" should be updated by the proxy
|
|
it = it_restart = 0;
|
|
it_restart_from_state_file = true;
|
|
|
|
// open the configfile
|
|
config_s.open (config_filename);
|
|
if (!config_s)
|
|
cvm::fatal_error ("Error: in opening configuration file \""+
|
|
std::string (config_filename)+"\".\n");
|
|
|
|
// read the config file into a string
|
|
std::string conf = "";
|
|
{
|
|
std::string line;
|
|
while (colvarparse::getline_nocomments (config_s, line))
|
|
conf.append (line+"\n");
|
|
// don't need the stream any more
|
|
config_s.close();
|
|
}
|
|
|
|
parse->get_keyval (conf, "analysis", b_analysis, false);
|
|
|
|
parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
|
|
colvarparse::parse_silent);
|
|
|
|
parse->get_keyval (conf, "eigenvalueCrossingThreshold",
|
|
colvarmodule::rotation::crossing_threshold, 1.0e-04,
|
|
colvarparse::parse_silent);
|
|
|
|
parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100);
|
|
parse->get_keyval (conf, "colvarsRestartFrequency", restart_out_freq,
|
|
proxy->restart_frequency());
|
|
|
|
// by default overwrite the existing trajectory file
|
|
parse->get_keyval (conf, "colvarsTrajAppend", cv_traj_append, false);
|
|
|
|
// input restart file
|
|
restart_in_name = proxy->input_prefix().size() ?
|
|
std::string (proxy->input_prefix()+".colvars.state") :
|
|
std::string ("") ;
|
|
|
|
// output restart file
|
|
restart_out_name = proxy->restart_output_prefix().size() ?
|
|
std::string (proxy->restart_output_prefix()+".colvars.state") :
|
|
std::string ("");
|
|
|
|
if (restart_out_name.size())
|
|
cvm::log ("The restart output state file will be \""+restart_out_name+"\".\n");
|
|
|
|
output_prefix = proxy->output_prefix();
|
|
|
|
cvm::log ("The final output state file will be \""+
|
|
(output_prefix.size() ?
|
|
std::string (output_prefix+".colvars.state") :
|
|
std::string ("colvars.state"))+"\".\n");
|
|
|
|
cv_traj_name =
|
|
(output_prefix.size() ?
|
|
std::string (output_prefix+".colvars.traj") :
|
|
std::string ("colvars.traj"));
|
|
cvm::log ("The trajectory file will be \""+
|
|
cv_traj_name+"\".\n");
|
|
|
|
if (cv_traj_freq) {
|
|
// open trajectory file
|
|
if (cv_traj_append) {
|
|
cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
|
|
"\".\n");
|
|
cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
|
|
} else {
|
|
proxy->backup_file (cv_traj_name.c_str());
|
|
cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
|
|
}
|
|
cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
|
|
}
|
|
|
|
// parse the options for collective variables
|
|
init_colvars (conf);
|
|
|
|
// parse the options for biases
|
|
init_biases (conf);
|
|
|
|
// done with the parsing, check that all keywords are valid
|
|
parse->check_keywords (conf, "colvarmodule");
|
|
cvm::log (cvm::line_marker);
|
|
|
|
// read the restart configuration, if available
|
|
if (restart_in_name.size()) {
|
|
// read the restart file
|
|
std::ifstream input_is (restart_in_name.c_str());
|
|
if (!input_is.good())
|
|
fatal_error ("Error: in opening restart file \""+
|
|
std::string (restart_in_name)+"\".\n");
|
|
else {
|
|
cvm::log ("Restarting from file \""+restart_in_name+"\".\n");
|
|
read_restart (input_is);
|
|
cvm::log (cvm::line_marker);
|
|
}
|
|
}
|
|
|
|
// check if it is possible to save output configuration
|
|
if ((!output_prefix.size()) && (!restart_out_name.size())) {
|
|
cvm::fatal_error ("Error: neither the final output state file or "
|
|
"the output restart file could be defined, exiting.\n");
|
|
}
|
|
|
|
cvm::log ("Collective variables module initialized.\n");
|
|
cvm::log (cvm::line_marker);
|
|
}
|
|
|
|
|
|
std::istream & colvarmodule::read_restart (std::istream &is)
|
|
{
|
|
{
|
|
// read global restart information
|
|
std::string restart_conf;
|
|
if (is >> colvarparse::read_block ("configuration", restart_conf)) {
|
|
if (it_restart_from_state_file) {
|
|
parse->get_keyval (restart_conf, "step",
|
|
it_restart, (size_t) 0,
|
|
colvarparse::parse_silent);
|
|
it = it_restart;
|
|
}
|
|
}
|
|
is.clear();
|
|
}
|
|
|
|
// colvars restart
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
if ( !((*cvi)->read_restart (is)) )
|
|
cvm::fatal_error ("Error: in reading restart configuration for collective variable \""+
|
|
(*cvi)->name+"\".\n");
|
|
}
|
|
|
|
// biases restart
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
if (!((*bi)->read_restart (is)))
|
|
fatal_error ("Error: in reading restart configuration for bias \""+
|
|
(*bi)->name+"\".\n");
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
return is;
|
|
}
|
|
|
|
|
|
|
|
void colvarmodule::init_colvars (std::string const &conf)
|
|
{
|
|
if (cvm::debug())
|
|
cvm::log ("Initializing the collective variables.\n");
|
|
|
|
std::string colvar_conf = "";
|
|
size_t pos = 0;
|
|
while (parse->key_lookup (conf, "colvar", colvar_conf, pos)) {
|
|
|
|
if (colvar_conf.size()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::increase_depth();
|
|
colvars.push_back (new colvar (colvar_conf));
|
|
(colvars.back())->check_keywords (colvar_conf, "colvar");
|
|
cvm::decrease_depth();
|
|
} else {
|
|
cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n");
|
|
}
|
|
colvar_conf = "";
|
|
}
|
|
|
|
|
|
if (!colvars.size())
|
|
cvm::fatal_error ("Error: no collective variables defined.\n");
|
|
|
|
if (colvars.size())
|
|
cvm::log (cvm::line_marker);
|
|
cvm::log ("Collective variables initialized, "+
|
|
cvm::to_str (colvars.size())+
|
|
" in total.\n");
|
|
}
|
|
|
|
|
|
void colvarmodule::init_biases (std::string const &conf)
|
|
{
|
|
if (cvm::debug())
|
|
cvm::log ("Initializing the collective variables biases.\n");
|
|
|
|
{
|
|
/// initialize ABF instances
|
|
std::string abf_conf = "";
|
|
size_t abf_pos = 0;
|
|
while (parse->key_lookup (conf, "abf", abf_conf, abf_pos)) {
|
|
if (abf_conf.size()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::increase_depth();
|
|
biases.push_back (new colvarbias_abf (abf_conf, "abf"));
|
|
(biases.back())->check_keywords (abf_conf, "abf");
|
|
cvm::decrease_depth();
|
|
n_abf_biases++;
|
|
} else {
|
|
cvm::log ("Warning: \"abf\" keyword found without configuration.\n");
|
|
}
|
|
abf_conf = "";
|
|
}
|
|
}
|
|
|
|
{
|
|
/// initialize harmonic restraints
|
|
std::string harm_conf = "";
|
|
size_t harm_pos = 0;
|
|
while (parse->key_lookup (conf, "harmonic", harm_conf, harm_pos)) {
|
|
if (harm_conf.size()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::increase_depth();
|
|
biases.push_back (new colvarbias_harmonic (harm_conf, "harmonic"));
|
|
(biases.back())->check_keywords (harm_conf, "harmonic");
|
|
cvm::decrease_depth();
|
|
n_harm_biases++;
|
|
} else {
|
|
cvm::log ("Warning: \"harmonic\" keyword found without configuration.\n");
|
|
}
|
|
harm_conf = "";
|
|
}
|
|
}
|
|
|
|
{
|
|
/// initialize histograms
|
|
std::string histo_conf = "";
|
|
size_t histo_pos = 0;
|
|
while (parse->key_lookup (conf, "histogram", histo_conf, histo_pos)) {
|
|
if (histo_conf.size()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::increase_depth();
|
|
biases.push_back (new colvarbias_histogram (histo_conf, "histogram"));
|
|
(biases.back())->check_keywords (histo_conf, "histogram");
|
|
cvm::decrease_depth();
|
|
n_histo_biases++;
|
|
} else {
|
|
cvm::log ("Warning: \"histogram\" keyword found without configuration.\n");
|
|
}
|
|
histo_conf = "";
|
|
}
|
|
}
|
|
|
|
{
|
|
/// initialize metadynamics instances
|
|
std::string meta_conf = "";
|
|
size_t meta_pos = 0;
|
|
while (parse->key_lookup (conf, "metadynamics", meta_conf, meta_pos)) {
|
|
if (meta_conf.size()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::increase_depth();
|
|
biases.push_back (new colvarbias_meta (meta_conf, "metadynamics"));
|
|
(biases.back())->check_keywords (meta_conf, "metadynamics");
|
|
cvm::decrease_depth();
|
|
n_meta_biases++;
|
|
} else {
|
|
cvm::log ("Warning: \"metadynamics\" keyword found without configuration.\n");
|
|
}
|
|
meta_conf = "";
|
|
}
|
|
}
|
|
|
|
if (biases.size())
|
|
cvm::log (cvm::line_marker);
|
|
cvm::log ("Collective variables biases initialized, "+
|
|
cvm::to_str (biases.size())+" in total.\n");
|
|
}
|
|
|
|
|
|
void colvarmodule::change_configuration(
|
|
std::string const &name, std::string const &conf)
|
|
{
|
|
cvm::increase_depth();
|
|
int found = 0;
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
if ( (*bi)->name == name ) {
|
|
++found;
|
|
(*bi)->change_configuration(conf);
|
|
}
|
|
}
|
|
if ( found < 1 ) cvm::fatal_error ("Error: bias not found");
|
|
if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name");
|
|
cvm::decrease_depth();
|
|
}
|
|
|
|
cvm::real colvarmodule::energy_difference(
|
|
std::string const &name, std::string const &conf)
|
|
{
|
|
cvm::increase_depth();
|
|
cvm::real energy_diff = 0.;
|
|
int found = 0;
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
if ( (*bi)->name == name ) {
|
|
++found;
|
|
energy_diff = (*bi)->energy_difference(conf);
|
|
}
|
|
}
|
|
if ( found < 1 ) cvm::fatal_error ("Error: bias not found");
|
|
if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name");
|
|
cvm::decrease_depth();
|
|
return energy_diff;
|
|
}
|
|
|
|
|
|
void colvarmodule::calc() {
|
|
cvm::real total_bias_energy = 0.0;
|
|
cvm::real total_colvar_energy = 0.0;
|
|
|
|
if (cvm::debug()) {
|
|
cvm::log (cvm::line_marker);
|
|
cvm::log ("Collective variables module, step no. "+
|
|
cvm::to_str (cvm::step_absolute())+"\n");
|
|
}
|
|
|
|
// calculate collective variables and their gradients
|
|
if (cvm::debug())
|
|
cvm::log ("Calculating collective variables.\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->calc();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
// update the biases and communicate their forces to the collective
|
|
// variables
|
|
if (cvm::debug() && biases.size())
|
|
cvm::log ("Updating collective variable biases.\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
total_bias_energy += (*bi)->update();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
// sum the forces from all biases for each collective variable
|
|
if (cvm::debug() && biases.size())
|
|
cvm::log ("Collecting forces from all biases.\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->reset_bias_force();
|
|
}
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
(*bi)->communicate_forces();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
if (cvm::b_analysis) {
|
|
// perform runtime analysis of colvars and biases
|
|
if (cvm::debug() && biases.size())
|
|
cvm::log ("Perform runtime analyses.\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->analyse();
|
|
}
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
(*bi)->analyse();
|
|
}
|
|
cvm::decrease_depth();
|
|
}
|
|
|
|
// sum up the forces for each colvar and integrate any internal
|
|
// equation of motion
|
|
if (cvm::debug())
|
|
cvm::log ("Updating the internal degrees of freedom "
|
|
"of colvars (if they have any).\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
total_colvar_energy += (*cvi)->update();
|
|
}
|
|
cvm::decrease_depth();
|
|
proxy->add_energy (total_bias_energy + total_colvar_energy);
|
|
|
|
// make collective variables communicate their forces to their
|
|
// coupled degrees of freedom (i.e. atoms)
|
|
if (cvm::debug())
|
|
cvm::log ("Communicating forces from the colvars to the atoms.\n");
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
if ((*cvi)->tasks[colvar::task_gradients])
|
|
(*cvi)->communicate_forces();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
// write restart file, if needed
|
|
if (restart_out_freq && !cvm::b_analysis) {
|
|
if ( (cvm::step_relative() > 0) &&
|
|
((cvm::step_absolute() % restart_out_freq) == 0) ) {
|
|
cvm::log ("Writing the state file \""+
|
|
restart_out_name+"\".\n");
|
|
proxy->backup_file (restart_out_name.c_str());
|
|
restart_out_os.open (restart_out_name.c_str());
|
|
restart_out_os.setf (std::ios::scientific, std::ios::floatfield);
|
|
if (!write_restart (restart_out_os))
|
|
cvm::fatal_error ("Error: in writing restart file.\n");
|
|
restart_out_os.close();
|
|
}
|
|
}
|
|
|
|
// write trajectory file, if needed
|
|
if (cv_traj_freq) {
|
|
|
|
if (cvm::debug())
|
|
cvm::log ("Writing trajectory file.\n");
|
|
|
|
// (re)open trajectory file
|
|
if (!cv_traj_os.good()) {
|
|
if (cv_traj_append) {
|
|
cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
|
|
"\".\n");
|
|
cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
|
|
} else {
|
|
cvm::log ("Overwriting colvar trajectory file \""+cv_traj_name+
|
|
"\".\n");
|
|
proxy->backup_file (cv_traj_name.c_str());
|
|
cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
|
|
}
|
|
cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
|
|
}
|
|
|
|
// write labels in the traj file every 1000 lines and at first ts
|
|
cvm::increase_depth();
|
|
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
|
|
cv_traj_os << "# " << cvm::wrap_string ("step", cvm::it_width-2)
|
|
<< " ";
|
|
if (cvm::debug())
|
|
cv_traj_os.flush();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->write_traj_label (cv_traj_os);
|
|
}
|
|
cv_traj_os << "\n";
|
|
if (cvm::debug())
|
|
cv_traj_os.flush();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
// write collective variable values to the traj file
|
|
cvm::increase_depth();
|
|
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
|
|
cv_traj_os << std::setw (cvm::it_width) << it
|
|
<< " ";
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->write_traj (cv_traj_os);
|
|
}
|
|
cv_traj_os << "\n";
|
|
if (cvm::debug())
|
|
cv_traj_os.flush();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
if (restart_out_freq) {
|
|
// flush the trajectory file if we are at the restart frequency
|
|
if ( (cvm::step_relative() > 0) &&
|
|
((cvm::step_absolute() % restart_out_freq) == 0) ) {
|
|
cvm::log ("Synchronizing (emptying the buffer of) trajectory file \""+
|
|
cv_traj_name+"\".\n");
|
|
cv_traj_os.flush();
|
|
}
|
|
}
|
|
} // end if (cv_traj_freq)
|
|
}
|
|
|
|
|
|
void colvarmodule::analyze()
|
|
{
|
|
if (cvm::debug()) {
|
|
cvm::log ("colvarmodule::analyze(), step = "+cvm::to_str (it)+".\n");
|
|
}
|
|
|
|
if (cvm::step_relative() == 0)
|
|
cvm::log ("Performing analysis.\n");
|
|
|
|
// perform colvar-specific analysis
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
cvm::increase_depth();
|
|
(*cvi)->analyse();
|
|
cvm::decrease_depth();
|
|
}
|
|
|
|
// perform bias-specific analysis
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
cvm::increase_depth();
|
|
(*bi)->analyse();
|
|
cvm::decrease_depth();
|
|
}
|
|
|
|
}
|
|
|
|
|
|
colvarmodule::~colvarmodule()
|
|
{
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
delete *cvi;
|
|
}
|
|
colvars.clear();
|
|
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
delete *bi;
|
|
}
|
|
biases.clear();
|
|
|
|
if (cv_traj_os.good()) {
|
|
cv_traj_os.close();
|
|
}
|
|
|
|
delete parse;
|
|
proxy = NULL;
|
|
}
|
|
|
|
|
|
void colvarmodule::write_output_files()
|
|
{
|
|
// if this is a simulation run (i.e. not a postprocessing), output data
|
|
// must be written to be able to restart the simulation
|
|
std::string const out_name =
|
|
(output_prefix.size() ?
|
|
std::string (output_prefix+".colvars.state") :
|
|
std::string ("colvars.state"));
|
|
cvm::log ("Saving collective variables state to \""+out_name+"\".\n");
|
|
proxy->backup_file (out_name.c_str());
|
|
std::ofstream out (out_name.c_str());
|
|
out.setf (std::ios::scientific, std::ios::floatfield);
|
|
this->write_restart (out);
|
|
out.close();
|
|
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->write_output_files();
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
// do not close to avoid problems with multiple NAMD runs
|
|
cv_traj_os.flush();
|
|
}
|
|
|
|
|
|
|
|
bool colvarmodule::read_traj (char const *traj_filename,
|
|
size_t traj_read_begin,
|
|
size_t traj_read_end)
|
|
{
|
|
cvm::log ("Opening trajectory file \""+
|
|
std::string (traj_filename)+"\".\n");
|
|
std::ifstream traj_is (traj_filename);
|
|
|
|
while (true) {
|
|
while (true) {
|
|
|
|
std::string line ("");
|
|
|
|
do {
|
|
if (!colvarparse::getline_nocomments (traj_is, line)) {
|
|
cvm::log ("End of file \""+std::string (traj_filename)+
|
|
"\" reached, or corrupted file.\n");
|
|
traj_is.close();
|
|
return false;
|
|
}
|
|
} while (line.find_first_not_of (colvarparse::white_space) == std::string::npos);
|
|
|
|
std::istringstream is (line);
|
|
|
|
if (!(is >> it)) return false;
|
|
|
|
if ( (it < traj_read_begin) ) {
|
|
|
|
if ((it % 1000) == 0)
|
|
std::cerr << "Skipping trajectory step " << it
|
|
<< " \r";
|
|
|
|
continue;
|
|
|
|
} else {
|
|
|
|
if ((it % 1000) == 0)
|
|
std::cerr << "Reading from trajectory, step = " << it
|
|
<< " \r";
|
|
|
|
if ( (traj_read_end > traj_read_begin) &&
|
|
(it > traj_read_end) ) {
|
|
std::cerr << "\n";
|
|
cvm::log ("Reached the end of the trajectory, "
|
|
"read_end = "+cvm::to_str (traj_read_end)+"\n");
|
|
return false;
|
|
}
|
|
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
if (!(*cvi)->read_traj (is)) {
|
|
cvm::log ("Error: in reading colvar \""+(*cvi)->name+
|
|
"\" from trajectory file \""+
|
|
std::string (traj_filename)+"\".\n");
|
|
return false;
|
|
}
|
|
}
|
|
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
std::ostream & colvarmodule::write_restart (std::ostream &os)
|
|
{
|
|
os << "configuration {\n"
|
|
<< " step " << std::setw (it_width)
|
|
<< it << "\n"
|
|
<< " dt " << dt() << "\n"
|
|
<< "}\n\n";
|
|
|
|
cvm::increase_depth();
|
|
for (std::vector<colvar *>::iterator cvi = colvars.begin();
|
|
cvi != colvars.end();
|
|
cvi++) {
|
|
(*cvi)->write_restart (os);
|
|
}
|
|
|
|
for (std::vector<colvarbias *>::iterator bi = biases.begin();
|
|
bi != biases.end();
|
|
bi++) {
|
|
(*bi)->write_restart (os);
|
|
}
|
|
cvm::decrease_depth();
|
|
|
|
return os;
|
|
}
|
|
|
|
|
|
|
|
void cvm::log (std::string const &message)
|
|
{
|
|
if (depth > 0)
|
|
proxy->log ((std::string (2*depth, ' '))+message);
|
|
else
|
|
proxy->log (message);
|
|
}
|
|
|
|
void cvm::increase_depth()
|
|
{
|
|
depth++;
|
|
}
|
|
|
|
void cvm::decrease_depth()
|
|
{
|
|
if (depth) depth--;
|
|
}
|
|
|
|
void cvm::fatal_error (std::string const &message)
|
|
{
|
|
proxy->fatal_error (message);
|
|
}
|
|
|
|
void cvm::exit (std::string const &message)
|
|
{
|
|
proxy->exit (message);
|
|
}
|
|
|
|
|
|
|
|
// static pointers
|
|
std::vector<colvar *> colvarmodule::colvars;
|
|
std::vector<colvarbias *> colvarmodule::biases;
|
|
size_t colvarmodule::n_abf_biases = 0;
|
|
size_t colvarmodule::n_harm_biases = 0;
|
|
size_t colvarmodule::n_histo_biases = 0;
|
|
size_t colvarmodule::n_meta_biases = 0;
|
|
colvarproxy *colvarmodule::proxy = NULL;
|
|
|
|
|
|
// static runtime data
|
|
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03;
|
|
size_t colvarmodule::it = 0;
|
|
size_t colvarmodule::it_restart = 0;
|
|
size_t colvarmodule::restart_out_freq = 0;
|
|
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;
|
|
|
|
|
|
// file name prefixes
|
|
std::string colvarmodule::output_prefix = "";
|
|
std::string colvarmodule::input_prefix = "";
|
|
std::string colvarmodule::restart_in_name = "";
|
|
|
|
|
|
// i/o constants
|
|
size_t const colvarmodule::it_width = 12;
|
|
size_t const colvarmodule::cv_prec = 14;
|
|
size_t const colvarmodule::cv_width = 21;
|
|
size_t const colvarmodule::en_prec = 14;
|
|
size_t const colvarmodule::en_width = 21;
|
|
std::string const colvarmodule::line_marker =
|
|
"----------------------------------------------------------------------\n";
|
|
|
|
|
|
|
|
|
|
|
|
|
|
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. 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)
|
|
{
|
|
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, 4, S_eigval, S_eigvec, &jac_nrot);
|
|
eigsrt (S_eigval, S_eigvec, 4);
|
|
// jacobi saves eigenvectors by columns
|
|
transpose (S_eigvec, 4);
|
|
|
|
// normalize eigenvectors
|
|
for (size_t ie = 0; ie < 4; ie++) {
|
|
cvm::real norm2 = 0.0;
|
|
for (size_t i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2));
|
|
cvm::real const norm = std::sqrt (norm2);
|
|
for (size_t 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: 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
|
|
for (size_t 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 (size_t 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);
|
|
void jacobi(cvm::real **a, int n, 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::fatal_error ("Too many iterations in routine jacobi.\n");
|
|
}
|
|
#undef ROTATE
|
|
|
|
|
|
void eigsrt(cvm::real d[], cvm::real **v, int n)
|
|
{
|
|
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, int n)
|
|
{
|
|
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;
|
|
}
|
|
}
|
|
}
|
|
|