From e2899f5df937729f36842ca03e09f5032ea21331 Mon Sep 17 00:00:00 2001
From: Andrew White <white.d.andrew@gmail.com>
Date: Sun, 9 Apr 2017 22:12:27 -0400
Subject: [PATCH] EDS: Added ability to have other CVs/Actions as centers
---
.../eds/rt-eds-adaptive/eds_restart.reference | 1 -
regtest/eds/rt-eds-append-restart/eds_restart | 1 -
.../eds_restart.reference | 2 -
.../eds/rt-eds-covar/eds_restart.reference | 1 -
regtest/eds/rt-eds-dynamic-center/Makefile | 1 +
regtest/eds/rt-eds-dynamic-center/config | 5 +
.../eds_restart.reference | 10 ++
regtest/eds/rt-eds-dynamic-center/plumed.dat | 3 +
regtest/eds/rt-eds-ramp/eds_restart.reference | 1 -
src/eds/EDS.cpp | 91 ++++++++++++++-----
10 files changed, 87 insertions(+), 29 deletions(-)
create mode 100644 regtest/eds/rt-eds-dynamic-center/Makefile
create mode 100644 regtest/eds/rt-eds-dynamic-center/config
create mode 100644 regtest/eds/rt-eds-dynamic-center/eds_restart.reference
create mode 100644 regtest/eds/rt-eds-dynamic-center/plumed.dat
diff --git a/regtest/eds/rt-eds-adaptive/eds_restart.reference b/regtest/eds/rt-eds-adaptive/eds_restart.reference
index 897e86a5d..588f16333 100644
--- a/regtest/eds/rt-eds-adaptive/eds_restart.reference
+++ b/regtest/eds/rt-eds-adaptive/eds_restart.reference
@@ -3,7 +3,6 @@
#! SET update_period 1
#! SET seed 0
#! SET kbt 2.49433863
- 0 1 0 0 0 7.483015889999999 0.1496603178
0 1 0 0 0 7.483015889999999 0.1496603178
1 1 0 0 0 7.483015889999999 0.1496603178
2 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178
diff --git a/regtest/eds/rt-eds-append-restart/eds_restart b/regtest/eds/rt-eds-append-restart/eds_restart
index 897e86a5d..588f16333 100644
--- a/regtest/eds/rt-eds-append-restart/eds_restart
+++ b/regtest/eds/rt-eds-append-restart/eds_restart
@@ -3,7 +3,6 @@
#! SET update_period 1
#! SET seed 0
#! SET kbt 2.49433863
- 0 1 0 0 0 7.483015889999999 0.1496603178
0 1 0 0 0 7.483015889999999 0.1496603178
1 1 0 0 0 7.483015889999999 0.1496603178
2 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178
diff --git a/regtest/eds/rt-eds-append-restart/eds_restart.reference b/regtest/eds/rt-eds-append-restart/eds_restart.reference
index 25d7633a0..1bf951dad 100644
--- a/regtest/eds/rt-eds-append-restart/eds_restart.reference
+++ b/regtest/eds/rt-eds-append-restart/eds_restart.reference
@@ -3,7 +3,6 @@
#! SET update_period 1
#! SET seed 0
#! SET kbt 2.49433863
- 0 1 0 0 0 7.483015889999999 0.1496603178
0 1 0 0 0 7.483015889999999 0.1496603178
1 1 0 0 0 7.483015889999999 0.1496603178
2 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178
@@ -14,7 +13,6 @@
#! SET update_period 1
#! SET seed 0
#! SET kbt 2.49433863
- 0 1 -1.107575371820062 0 -1.107575371820062 7.483015889999999 0.1496603178
0 1 -1.107575371820062 0 -1.107575371820062 7.483015889999999 0.1496603178
1 1 -1.107575371820062 0 -1.107575371820062 7.483015889999999 0.1496603178
2 1 -2.215150743640125 0 -1.107575371820062 7.483015889999999 0.1496603178
diff --git a/regtest/eds/rt-eds-covar/eds_restart.reference b/regtest/eds/rt-eds-covar/eds_restart.reference
index 94b7bf8a2..0ff3608df 100644
--- a/regtest/eds/rt-eds-covar/eds_restart.reference
+++ b/regtest/eds/rt-eds-covar/eds_restart.reference
@@ -3,7 +3,6 @@
#! SET update_period 1
#! SET seed 574
#! SET kbt 2.49433863
- 0 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178
0 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178
1 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178 1 0 0 0 7.483015889999999 0.1496603178
2 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178 1 -1.107575371820062 0 0 7.483015889999999 0.1496603178
diff --git a/regtest/eds/rt-eds-dynamic-center/Makefile b/regtest/eds/rt-eds-dynamic-center/Makefile
new file mode 100644
index 000000000..3703b27ce
--- /dev/null
+++ b/regtest/eds/rt-eds-dynamic-center/Makefile
@@ -0,0 +1 @@
+include ../../scripts/test.make
diff --git a/regtest/eds/rt-eds-dynamic-center/config b/regtest/eds/rt-eds-dynamic-center/config
new file mode 100644
index 000000000..d524a9db4
--- /dev/null
+++ b/regtest/eds/rt-eds-dynamic-center/config
@@ -0,0 +1,5 @@
+type=driver
+plumed_modules=eds
+arg="--plumed plumed.dat --ixyz trajectory.xyz"
+extra_files="../../trajectories/trajectory.xyz"
+
diff --git a/regtest/eds/rt-eds-dynamic-center/eds_restart.reference b/regtest/eds/rt-eds-dynamic-center/eds_restart.reference
new file mode 100644
index 000000000..5d98f4ca4
--- /dev/null
+++ b/regtest/eds/rt-eds-dynamic-center/eds_restart.reference
@@ -0,0 +1,10 @@
+#! FIELDS time d1_center d1_set d1_target d1_coupling d1_maxrange d1_maxgrad
+#! SET adaptive 1
+#! SET update_period 1
+#! SET seed 0
+#! SET kbt 2.49433863
+ 0 1.097204685252765 0 0 0 7.483015889999999 0.1496603178
+ 1 1.058762371055742 0 0 0 7.483015889999999 0.1496603178
+ 2 1.095818674401286 -1.107575371820062 0 0 7.483015889999999 0.1496603178
+ 3 1.162848545384402 -1.107575371820062 0 -1.107575371820062 7.483015889999999 0.1496603178
+ 4 1.21603261267908 -1.107575371820062 0 -1.107575371820062 7.483015889999999 0.1496603178
diff --git a/regtest/eds/rt-eds-dynamic-center/plumed.dat b/regtest/eds/rt-eds-dynamic-center/plumed.dat
new file mode 100644
index 000000000..b1c732129
--- /dev/null
+++ b/regtest/eds/rt-eds-dynamic-center/plumed.dat
@@ -0,0 +1,3 @@
+d2: DISTANCE ATOMS=3,4
+d1: DISTANCE ATOMS=1,2
+EDS ARG=d1 CENTER_ARG=d2 BIAS_SCALE=1.0 PERIOD=2 LABEL=eds OUT_RESTART=eds_restart TEMP=300.0
diff --git a/regtest/eds/rt-eds-ramp/eds_restart.reference b/regtest/eds/rt-eds-ramp/eds_restart.reference
index b73cc9fba..15e9b1cd2 100644
--- a/regtest/eds/rt-eds-ramp/eds_restart.reference
+++ b/regtest/eds/rt-eds-ramp/eds_restart.reference
@@ -3,7 +3,6 @@
#! SET update_period -2
#! SET seed 0
#! SET kbt 2.49433863
- 0 1 0.5 -3 0.5 7.483015889999999 0.1496603178
0 1 0.5 -3 0.5 7.483015889999999 0.1496603178
1 1 0.5 -3 -1.25 7.483015889999999 0.1496603178
2 1 0.5 -3 -3 7.483015889999999 0.1496603178
diff --git a/src/eds/EDS.cpp b/src/eds/EDS.cpp
index 0ba032c95..6d80bca5f 100644
--- a/src/eds/EDS.cpp
+++ b/src/eds/EDS.cpp
@@ -111,6 +111,7 @@ private:
/*We will get this and store it once, since on-the-fly changing number of CVs will be fatal*/
const unsigned int ncvs_;
std::vector<double> center_;
+ std::vector<Value*> center_values_;
std::vector<double> scale_;
std::vector<double> current_coupling_;
std::vector<double> set_coupling_;
@@ -128,6 +129,7 @@ private:
std::string out_restart_name_;
OFile out_restart_;
IFile in_restart_;
+ bool b_c_values_;
bool b_adaptive_;
bool b_freeze_;
bool b_equil_;
@@ -173,7 +175,10 @@ public:
void EDS::registerKeywords(Keywords& keys){
Bias::registerKeywords(keys);
keys.use("ARG");
- keys.add("compulsory","CENTER","The desired centers (equilibrium values) which will be sought during the adaptive linear biasing.");
+ keys.add("optional","CENTER","The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. This is for fixed values");
+ keys.add("optional","CENTER_ARG","The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. "
+ "CENTER_ARG is for calculated centers, e.g. from a CV or analysis. ");
+
keys.add("compulsory","PERIOD","Steps over which to adjust bias");
keys.add("compulsory","RANGE","3.0","The largest magnitude of the force constant which one expects (in kBT) for each CV based");
@@ -208,7 +213,6 @@ public:
EDS::EDS(const ActionOptions&ao):
PLUMED_BIAS_INIT(ao),
ncvs_(getNumberOfArguments()),
- center_(ncvs_,1.0),
scale_(ncvs_,0.0),
current_coupling_(ncvs_,0.0),
set_coupling_(ncvs_,0.0),
@@ -254,6 +258,7 @@ public:
}
parseVector("CENTER",center_);
+ parseArgumentList("CENTER_ARG",center_values_);
parseVector("BIAS_SCALE", scale_);
parseVector("RANGE",max_coupling_range_);
parseVector("FIXED",target_coupling_);
@@ -270,13 +275,44 @@ public:
parse("IN_RESTART",in_restart_name_);
checkRead();
+ /*
+ * Things that are different when usnig chaning centers:
+ * 1. Scale
+ * 2. The log file
+ * 3. Reading Restarts
+ */
+
+ if(center_.size() == 0) {
+ if(center_values_.size() == 0)
+ error("Must set either CENTER or CENTER_ARG");
+ else if(center_values_.size() != ncvs_)
+ error("CENTER_ARG must contain the same number of variables as ARG");
+ b_c_values_ = true;
+ center_.resize(ncvs_);
+ log.printf(" EDS will use possibly varying centers\n");
+ } else {
+ if(center_.size() != ncvs_)
+ error("Must have same number of CENTER arguments as ARG arguments");
+ else if(center_values_.size() != 0)
+ error("You can only set CENTER or CENTER_ARG. Not both");
+ b_c_values_ = false;
+ log.printf(" EDS will use fixed centers\n");
+ }
+
+
+
log.printf(" setting scaling:");
if(scale_.size() > 0 && scale_.size() < ncvs_) {
error("the number of BIAS_SCALE values be the same as number of CVs");
- } else if(scale_.size() == 0) {
+ } else if(scale_.size() == 0 && b_c_values_) {
+ log.printf(" Setting SCALE to be 1 for all CVs\n");
+ scale_.resize(ncvs_);
+ for(unsigned int i = 0; i < ncvs_; ++i)
+ scale_[i] = 1;
+ } else if(scale_.size() == 0 && !b_c_values_) {
log.printf(" (default) ");
- scale_.resize(center_.size());
+ scale_.resize(ncvs_);
for(unsigned int i = 0; i < scale_.size(); i++) {
if(center_[i]==0)
error("BIAS_SCALE parameter has been set to CENTER value of 0 (as is default). This will divide by 0, so giving up. See doc for EDS bias");
@@ -322,14 +358,20 @@ public:
log.printf(" kBT = %f\n",kbt_);
log.printf(" Updating every %i steps\n",update_period_);
- log.printf(" with centers:");
- for(unsigned int i = 0;i<center_.size();i++){
- log.printf(" %f",center_[i]);
+ if(!b_c_values_) {
+ log.printf(" with centers:");
+ for(unsigned int i = 0;i< ncvs_;i++){
+ log.printf(" %f ",center_[i]);
+ }
+ } else {
+ log.printf(" with actions centers:");
+ for(unsigned int i = 0;i< ncvs_;i++){
+ log.printf(" %s ",center_values_[i]->getName().c_str());
+ //add dependency on these actions
+ addDependency(center_values_[i]->getPntrToAction());
+ }
}
- for(unsigned int i = 0; i < scale_.size(); i++)
- log.printf(" %f",center_[i]);
-
log.printf("\n with initial ranges / rates:\n");
for(unsigned int i = 0;i<max_coupling_range_.size();i++) {
//this is just an empirical guess. Bigger range, bigger grads. Less frequent updates, bigger changes
@@ -448,7 +490,7 @@ public:
for(unsigned int i = 0;i<ncvs_;++i) {
cv_name = getPntrToArgument(i)->getName();
- in_restart_.scanField(cv_name + +"_center",center_[i]);
+ in_restart_.scanField(cv_name + "_center", set_coupling_[i]);
in_restart_.scanField(cv_name + "_set", set_coupling_[i]);
in_restart_.scanField(cv_name + "_target",target_coupling_[i]);
in_restart_.scanField(cv_name + "_coupling",current_coupling_[i]);
@@ -510,7 +552,6 @@ public:
out_restart_.addConstantField("seed").printField("seed",seed_);
out_restart_.addConstantField("kbt").printField("kbt",kbt_);
- writeOutRestart();
}
void EDS::writeOutRestart() {
@@ -519,7 +560,6 @@ public:
for(unsigned int i = 0;i<ncvs_;++i) {
cv_name = getPntrToArgument(i)->getName();
-
out_restart_.printField(cv_name + "_center",center_[i]);
out_restart_.printField(cv_name + "_set",set_coupling_[i]);
out_restart_.printField(cv_name + "_target",target_coupling_[i]);
@@ -533,6 +573,11 @@ public:
void EDS::calculate(){
+
+ //get center values from action if necessary
+ if(b_c_values_)
+ for(unsigned int i = 0; i < ncvs_; ++i)
+ center_[i] = center_values_[i]->get();
apply_bias();
@@ -542,10 +587,11 @@ public:
//check if we're ramping or doing normal updates and then restart if needed. The ramping check
//is complicated because we could be frozen, finished ramping or not ramping.
//The + 2 is so we have an extra line showing that the bias isn't changing (for my sanity and yours)
- if( b_write_restart_ &&
+ if( b_write_restart_){
+ if(getStep() == 0 ||
( (update_period_ < 0 && !b_freeze_ && update_calls_ <= fabs(update_period_) + 2) ||
- (update_period_ > 0 && update_calls_ % update_period_ == 0 ) ) ) {
- writeOutRestart();
+ (update_period_ > 0 && update_calls_ % update_period_ == 0 ) ) )
+ writeOutRestart();
}
int b_finished_equil_flag = 1;
@@ -613,13 +659,12 @@ public:
void EDS::apply_bias() {
//Compute linear force as in "restraint"
- double ene = 0;
- double totf2 = 0;
+ double ene = 0, totf2 = 0, cv, m , f;
for(unsigned int i = 0; i < ncvs_; ++i) {
- const double cv = difference(i, center_[i], getArgument(i));
- const double m = current_coupling_[i];
- const double f = -m;
+ cv = difference(i, center_[i], getArgument(i));
+ m = current_coupling_[i];
+ f = -m;
ene += m*cv;
setOutputForce(i,f);
totf2 += f*f;
@@ -671,7 +716,7 @@ public:
for(unsigned int i = 0; i< ncvs_; ++i){
tmp = 0;
for(unsigned int j = 0; j < ncvs_; ++j)
- tmp += (means_[i] - center_[i]) * covar_(i,j);
+ tmp += difference(i, means_[i], center_[i]) * covar_(i,j);
step_size_[i] = 2 * tmp / kbt_ / scale_[i] * update_calls_ / (update_calls_ - 1);
}
@@ -680,7 +725,7 @@ public:
void EDS::calc_ssd_step_size() {
double tmp;
for(unsigned int i = 0; i< ncvs_; ++i){
- tmp = 2. * (means_[i] - center_[i]) * ssds_[i] / (update_calls_ - 1);
+ tmp = 2. * difference(i, means_[i], center_[i]) * ssds_[i] / (update_calls_ - 1);
step_size_[i] = tmp / kbt_/scale_[i];
}
}
--
GitLab