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

update colvars library to version 2016-08-05 

(cherry picked from commit 459db2eb6b)

# Conflicts:
#	doc/src/PDF/colvars-refman-lammps.pdf
This commit is contained in:
Axel Kohlmeyer 2016-09-06 21:26:50 -04:00
parent 0252347d43
commit 77620106a4
7 changed files with 152 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
lib/colvars/colvar.cpp
View File

@ -1187,40 +1187,6 @@ cvm::real colvar::update_forces_energy()
f -= fj;
}
if (is_enabled(f_cv_extended_Lagrangian)) {
cvm::real dt = cvm::dt();
cvm::real f_ext;
// the total force is applied to the fictitious mass, while the
// atoms only feel the harmonic force
// fr: bias force on extended coordinate (without harmonic spring), for output in trajectory
// f_ext: total force on extended coordinate (including harmonic spring)
// f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
// (note: wall potential is added to f after this block)
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
// Because of leapfrog, kinetic energy at time i is approximate
kinetic_energy = 0.5 * ext_mass * vr * vr;
potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
// leap to v_(i+1/2)
if (is_enabled(f_cv_Langevin)) {
vr -= dt * ext_gamma * vr.real_value;
vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
}
vr += (0.5 * dt) * f_ext / ext_mass;
xr += dt * vr;
xr.apply_constraints();
if (this->b_periodic) this->wrap(xr);
}
// Adding wall potential to "true" colvar force, whether or not an extended coordinate is in use
if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) {
// Wall force
@ -1257,6 +1223,37 @@ cvm::real colvar::update_forces_energy()
}
}
if (is_enabled(f_cv_extended_Lagrangian)) {
cvm::real dt = cvm::dt();
cvm::real f_ext;
// the total force is applied to the fictitious mass, while the
// atoms only feel the harmonic force
// fr: bias force on extended coordinate (without harmonic spring), for output in trajectory
// f_ext: total force on extended coordinate (including harmonic spring)
// f: - initially, external biasing force
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
// Because of leapfrog, kinetic energy at time i is approximate
kinetic_energy = 0.5 * ext_mass * vr * vr;
potential_energy = 0.5 * ext_force_k * this->dist2(xr, x);
// leap to v_(i+1/2)
if (is_enabled(f_cv_Langevin)) {
vr -= dt * ext_gamma * vr.real_value;
vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
}
vr += (0.5 * dt) * f_ext / ext_mass;
xr += dt * vr;
xr.apply_constraints();
if (this->b_periodic) this->wrap(xr);
}
f_accumulated += f;
if (is_enabled(f_cv_fdiff_velocity)) {

78
lib/colvars/colvarbias_abf.cpp
View File

@ -12,6 +12,7 @@ colvarbias_abf::colvarbias_abf(char const *key)
samples(NULL),
z_gradients(NULL),
z_samples(NULL),
czar_gradients(NULL),
last_gradients(NULL),
last_samples(NULL)
{
@ -144,6 +145,7 @@ int colvarbias_abf::init(std::string const &conf)
z_gradients->request_actual_value();
z_gradients->samples = z_samples;
z_samples->has_parent_data = true;
czar_gradients = new colvar_grid_gradient(colvars);
}
// For shared ABF, we store a second set of grids.
@ -178,6 +180,21 @@ colvarbias_abf::~colvarbias_abf()
gradients = NULL;
}
if (z_samples) {
delete z_samples;
z_samples = NULL;
}
if (z_gradients) {
delete z_gradients;
z_gradients = NULL;
}
if (czar_gradients) {
delete czar_gradients;
czar_gradients = NULL;
}
// shared ABF
// We used to only do this if "shared" was defined,
// but now we can call shared externally
@ -237,7 +254,11 @@ int colvarbias_abf::update()
z_bin[i] = z_samples->current_bin_scalar(i);
}
if ( z_samples->index_ok(z_bin) ) {
// Set increment flag to 0 to only increment
for (size_t i=0; i<colvars.size(); i++) {
// If we are outside the range of xi, the force has not been obtained above
// the function is just an accessor, so cheap to call again anyway
force[i] = colvars[i]->system_force();
}
z_gradients->acc_force(z_bin, force);
}
}
@ -404,14 +425,10 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
{
std::string samples_out_name = prefix + ".count";
std::string gradients_out_name = prefix + ".grad";
std::string z_gradients_out_name = prefix + ".zgrad";
std::string z_samples_out_name = prefix + ".zcount";
std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
cvm::ofstream samples_os;
cvm::ofstream gradients_os;
cvm::ofstream z_samples_os;
cvm::ofstream z_gradients_os;
if (!append) cvm::backup_file(samples_out_name.c_str());
samples_os.open(samples_out_name.c_str(), mode);
@ -442,6 +459,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
}
if (z_gradients) {
// Write eABF-related quantities
std::string z_samples_out_name = prefix + ".zcount";
std::string z_gradients_out_name = prefix + ".zgrad";
std::string czar_gradients_out_name = prefix + ".czar";
cvm::ofstream z_samples_os;
cvm::ofstream z_gradients_os;
cvm::ofstream czar_gradients_os;
if (!append) cvm::backup_file(z_samples_out_name.c_str());
z_samples_os.open(z_samples_out_name.c_str(), mode);
if (!z_samples_os.is_open()) {
@ -458,6 +483,24 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
z_gradients->write_multicol(z_gradients_os);
z_gradients_os.close();
// Calculate CZAR estimator of gradients
for (std::vector<int> ix = czar_gradients->new_index();
czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
- cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n),
n);
}
}
if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
if (!czar_gradients_os.is_open()) {
cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
}
czar_gradients->write_multicol(czar_gradients_os);
czar_gradients_os.close();
if (colvars.size() == 1) {
std::string z_pmf_out_name = prefix + ".zpmf";
if (!append) cvm::backup_file(z_pmf_out_name.c_str());
@ -468,6 +511,16 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
z_gradients->write_1D_integral(z_pmf_os);
z_pmf_os << std::endl;
z_pmf_os.close();
std::string czar_pmf_out_name = prefix + ".czarpmf";
if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
cvm::ofstream czar_pmf_os;
// Do numerical integration and output a PMF
czar_pmf_os.open(czar_pmf_out_name.c_str(), mode);
if (!czar_pmf_os.is_open()) cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
czar_gradients->write_1D_integral(czar_pmf_os);
czar_pmf_os << std::endl;
czar_pmf_os.close();
}
}
@ -545,24 +598,27 @@ std::ostream & colvarbias_abf::write_restart(std::ostream& os)
{
std::ios::fmtflags flags(os.flags());
os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
os << "abf {\n"
<< " configuration {\n"
<< " name " << this->name << "\n";
os << " }\n";
os << "samples\n";
os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
os << "\nsamples\n";
samples->write_raw(os, 8);
os.flags(flags);
os << "\ngradient\n";
gradients->write_raw(os);
gradients->write_raw(os, 8);
if (z_gradients) {
os << "z_samples\n";
os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
os << "\nz_samples\n";
z_samples->write_raw(os, 8);
os.flags(flags);
os << "\nz_gradient\n";
z_gradients->write_raw(os);
z_gradients->write_raw(os, 8);
}
os << "}\n\n";
@ -638,7 +694,7 @@ std::istream & colvarbias_abf::read_restart(std::istream& is)
}
if (z_gradients) {
if ( !(is >> key) || !(key == "z_samples")) {
if ( !(is >> key) || !(key == "z_samples")) {
cvm::log("Error: in reading restart configuration for ABF bias \""+
this->name+"\" at position "+
cvm::to_str(is.tellg())+" in stream.\n");

2
lib/colvars/colvarbias_abf.h
View File

@ -60,6 +60,8 @@ private:
colvar_grid_gradient *z_gradients;
/// n-dim grid of number of samples on "real" coordinate for eABF z-based estimator
colvar_grid_count *z_samples;
/// n-dim grid contining CZAR estimator of "real" free energy gradients
colvar_grid_gradient *czar_gradients;
// shared ABF
bool shared_on;

7
lib/colvars/colvarbias_histogram.cpp
View File

@ -75,13 +75,12 @@ int colvarbias_histogram::init(std::string const &conf)
}
grid = new colvar_grid_scalar();
grid->init_from_colvars(colvars);
{
std::string grid_conf;
if (key_lookup(conf, "grid", grid_conf)) {
grid->parse_params(grid_conf);
} else {
grid->init_from_colvars(colvars);
}
}
@ -256,6 +255,9 @@ std::istream & colvarbias_histogram::read_restart(std::istream& is)
std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
{
std::ios::fmtflags flags(os.flags());
os.setf(std::ios::fmtflags(0), std::ios::floatfield);
os << "histogram {\n"
<< " configuration {\n"
<< " name " << this->name << "\n";
@ -266,5 +268,6 @@ std::ostream & colvarbias_histogram::write_restart(std::ostream& os)
os << "}\n\n";
os.flags(flags);
return os;
}

4
lib/colvars/colvargrid.cpp
View File

@ -35,13 +35,13 @@ colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
: colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
: colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL), grad(NULL)
{
grad = new cvm::real[nd];
}
colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
: colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
: colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL), grad(NULL)
{
grad = new cvm::real[nd];
}

44
lib/colvars/colvargrid.h
View File

@ -1160,6 +1160,50 @@ public:
}
has_data = true;
}
/// \brief Return the log-gradient from finite differences
/// on the *same* grid for dimension n
inline const cvm::real log_gradient_finite_diff( const std::vector<int> &ix0,
int n = 0)
{
cvm::real A0, A1;
std::vector<int> ix;
// factor for mesh width, 2.0 for central finite difference
// but only 1.0 on edges for non-PBC coordinates
cvm::real factor;
if (periodic[n]) {
factor = 2.;
ix = ix0;
ix[n]--; wrap(ix);
A0 = data[address(ix)];
ix = ix0;
ix[n]++; wrap(ix);
A1 = data[address(ix)];
} else {
factor = 0.;
ix = ix0;
if (ix[n] > 0) { // not left edge
ix[n]--;
factor += 1.;
}
A0 = data[address(ix)];
ix = ix0;
if (ix[n]+1 < nx[n]) { // not right edge
ix[n]++;
factor += 1.;
}
A1 = data[address(ix)];
}
if (A0 == 0 || A1 == 0) {
// can't handle empty bins
return 0.;
} else {
return (std::log((cvm::real)A1) - std::log((cvm::real)A0))
/ (widths[n] * factor);
}
}
};

2
lib/colvars/colvarmodule.h
View File

@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2016-08-01"
#define COLVARS_VERSION "2016-08-05"
#endif
#ifndef COLVARS_DEBUG