diff --git a/CHANGES/Unreleased.txt b/CHANGES/Unreleased.txt
index 86b7bb8201a2557c7cc01eb1bb973e26e5e41825..dc96c96c73cb52ae61973bc4d00c1a51131ce399 100644
--- a/CHANGES/Unreleased.txt
+++ b/CHANGES/Unreleased.txt
@@ -46,5 +46,6 @@ Changes from version 2.2 which are relevant for users:
- \ref COMMITTOR can now be used to define multiple basins
- The number of atoms admitted in \ref BRIDGE has been significantly increased, see \issue{185}.
- \ref driver now allows --trajectory-stride to be set to zero when reading with --ixtc/--itrr. In this case, step number is read from the trajectory file.
+ - MetaD and PBMetaD can now be restarted from a GRID
*/
diff --git a/CHANGES/v2.2.txt b/CHANGES/v2.2.txt
index 076495368c4a80d517f8bb6669bc3b4d532b470e..5d18c6bb4f0b9a6929510133737fe6faa74445e3 100644
--- a/CHANGES/v2.2.txt
+++ b/CHANGES/v2.2.txt
@@ -137,7 +137,7 @@ For developers:
in cmd strings.
- ./configure is not automatically relaunched anymore when doing `make clean`.
-Version 2.2.3 ()
+Unreleased changes (will be included in 2.2.3)
----------------------------------------------
For users:
diff --git a/regtest/basic/rt-mpi7/COLVAR.0.reference b/regtest/basic/rt-mpi7/COLVAR.0.reference
new file mode 100644
index 0000000000000000000000000000000000000000..1a5718c4dbbe04d43b649d74bab59466dd83f09f
--- /dev/null
+++ b/regtest/basic/rt-mpi7/COLVAR.0.reference
@@ -0,0 +1,11 @@
+#! FIELDS time d1 d2 meta1.bias meta2.bias restart.bias
+ 0.000000 1.262593 1.097205 0.000000 0.000000 0.810970
+ 0.050000 1.317587 1.058762 0.000000 0.000000 0.886420
+ 0.100000 1.393388 1.095819 0.000000 0.000000 1.059839
+ 0.150000 1.475479 1.162849 0.321330 0.321330 0.867021
+ 0.200000 1.490756 1.216033 0.248441 0.248441 0.701164
+ 0.250000 1.262593 1.097205 0.419631 0.419631 0.810970
+ 0.300000 1.317587 1.058762 0.460037 0.460037 0.886420
+ 0.350000 1.393388 1.095819 0.781166 0.781166 1.059839
+ 0.400000 1.475479 1.162849 0.593033 0.593033 0.867021
+ 0.450000 1.490756 1.216033 0.701164 0.701164 0.701164
diff --git a/regtest/basic/rt-mpi7/COLVAR.1.reference b/regtest/basic/rt-mpi7/COLVAR.1.reference
new file mode 100644
index 0000000000000000000000000000000000000000..05239af624b9f0c5c012c6e545e89777a812b2ae
--- /dev/null
+++ b/regtest/basic/rt-mpi7/COLVAR.1.reference
@@ -0,0 +1,11 @@
+#! FIELDS time d1 d2 meta1.bias meta2.bias restart.bias
+ 0.000000 1.490756 1.216033 0.000000 0.000000 0.608488
+ 0.050000 1.475479 1.162849 0.000000 0.000000 0.777038
+ 0.100000 1.393388 1.095819 0.000000 0.000000 1.015893
+ 0.150000 1.317587 1.058762 0.371518 0.371518 0.884771
+ 0.200000 1.262593 1.097205 0.333056 0.333056 0.817604
+ 0.250000 1.490756 1.216033 0.277303 0.277303 0.608488
+ 0.300000 1.475479 1.162849 0.363420 0.363420 0.777038
+ 0.350000 1.393388 1.095819 0.785000 0.785000 1.015893
+ 0.400000 1.317587 1.058762 0.665250 0.665250 0.884771
+ 0.450000 1.262593 1.097205 0.817604 0.817604 0.817604
diff --git a/regtest/basic/rt-mpi7/COLVAR.2.reference b/regtest/basic/rt-mpi7/COLVAR.2.reference
new file mode 100644
index 0000000000000000000000000000000000000000..7c9e1a2b57e65742301f519f23133a35ebf48517
--- /dev/null
+++ b/regtest/basic/rt-mpi7/COLVAR.2.reference
@@ -0,0 +1,11 @@
+#! FIELDS time d1 d2 meta1.bias meta2.bias restart.bias
+ 0.000000 1.317587 1.058762 0.000000 0.000000 0.748387
+ 0.050000 1.393388 1.095819 0.000000 0.000000 0.903920
+ 0.100000 1.262593 1.097205 0.000000 0.000000 0.696759
+ 0.150000 1.490756 1.216033 0.028881 0.028881 0.633495
+ 0.200000 1.475479 1.162849 0.042140 0.042140 0.769058
+ 0.250000 1.317587 1.058762 0.284901 0.284901 0.748387
+ 0.300000 1.393388 1.095819 0.347138 0.347138 0.903920
+ 0.350000 1.262593 1.097205 0.616792 0.616792 0.696759
+ 0.400000 1.490756 1.216033 0.512368 0.512368 0.633495
+ 0.450000 1.475479 1.162849 0.769058 0.769058 0.769058
diff --git a/regtest/basic/rt-mpi7/grid.0.reference b/regtest/basic/rt-mpi7/grid.0.reference
index eba63d858b83f3706f097b848ffd8a90a827381e..230fe5b5d9fd1b369d3369ef43304441be072aae 100644
--- a/regtest/basic/rt-mpi7/grid.0.reference
+++ b/regtest/basic/rt-mpi7/grid.0.reference
@@ -1,4 +1,4 @@
-#! FIELDS d1 d2 @3.bias der_d1 der_d2
+#! FIELDS d1 d2 meta1.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/grid.1.reference b/regtest/basic/rt-mpi7/grid.1.reference
index f68da103b2e18c7283d0d1de06dbb3659ce1fc01..5fe44963023632d4a3978d386b7c8b046fbe248b 100644
--- a/regtest/basic/rt-mpi7/grid.1.reference
+++ b/regtest/basic/rt-mpi7/grid.1.reference
@@ -1,4 +1,4 @@
-#! FIELDS d1 d2 @3.bias der_d1 der_d2
+#! FIELDS d1 d2 meta1.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/grid.2.reference b/regtest/basic/rt-mpi7/grid.2.reference
index eb0d399a96a6f8fb64d9b5836cbe5d707de4378e..840a36e466c9602667da1adbb964c3cb70b37b1c 100644
--- a/regtest/basic/rt-mpi7/grid.2.reference
+++ b/regtest/basic/rt-mpi7/grid.2.reference
@@ -1,4 +1,4 @@
-#! FIELDS d1 d2 @3.bias der_d1 der_d2
+#! FIELDS d1 d2 meta1.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/gridx.0.reference b/regtest/basic/rt-mpi7/gridx.0.reference
index 1dd64c8f88e986f942022283fe36f2ea32ae979f..624de6d66113cf0d3bd7134cb59cde19de0eab41 100644
--- a/regtest/basic/rt-mpi7/gridx.0.reference
+++ b/regtest/basic/rt-mpi7/gridx.0.reference
@@ -1,5 +1,5 @@
# this is to test restart
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -37,7 +37,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -75,7 +75,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -113,7 +113,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -151,7 +151,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/gridx.1.reference b/regtest/basic/rt-mpi7/gridx.1.reference
index f6f7385b4ce179cfac422857e1d998a0a8f24600..137826c6fe4fcfb3074f8ceb782ac7fc7bc63ee3 100644
--- a/regtest/basic/rt-mpi7/gridx.1.reference
+++ b/regtest/basic/rt-mpi7/gridx.1.reference
@@ -1,4 +1,4 @@
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -36,7 +36,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -74,7 +74,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -112,7 +112,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -150,7 +150,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/gridx.2.reference b/regtest/basic/rt-mpi7/gridx.2.reference
index aaccbbd58c3ac40362c64282c1f3644a834d13f1..9c8da6c0d6ff73e9bffb2e169300d57e5a4cf113 100644
--- a/regtest/basic/rt-mpi7/gridx.2.reference
+++ b/regtest/basic/rt-mpi7/gridx.2.reference
@@ -1,4 +1,4 @@
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -36,7 +36,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -74,7 +74,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -112,7 +112,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
@@ -150,7 +150,7 @@
1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
-#! FIELDS d1 d2 @4.bias der_d1 der_d2
+#! FIELDS d1 d2 meta2.bias der_d1 der_d2
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1 5
diff --git a/regtest/basic/rt-mpi7/plumed.dat b/regtest/basic/rt-mpi7/plumed.dat
index af5bbaacca43ed6bc758226f163a35fb849196a3..20044e1cab8d09878a8df1407316740fd1694f27 100644
--- a/regtest/basic/rt-mpi7/plumed.dat
+++ b/regtest/basic/rt-mpi7/plumed.dat
@@ -8,6 +8,7 @@ METAD ...
FMT=%6.2f
GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4
GRID_WFILE=grid GRID_WSTRIDE=2
+ LABEL=meta1
...
METAD ...
@@ -16,8 +17,16 @@ METAD ...
GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4
GRID_WFILE=gridx GRID_WSTRIDE=2 STORE_GRIDS
FILE=HILLSX
+ LABEL=meta2
...
+METAD ...
+ ARG=d1,d2 SIGMA=0.1,0.1 PACE=2 HEIGHT=0.0
+ FMT=%6.2f
+ GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4
+ GRID_RFILE=rgrid
+ FILE=HILLSR
+ LABEL=restart
+...
-
-
+PRINT ARG=d1,d2,meta1.bias,meta2.bias,restart.bias FILE=COLVAR
diff --git a/regtest/basic/rt-mpi7/rgrid.0 b/regtest/basic/rt-mpi7/rgrid.0
new file mode 100644
index 0000000000000000000000000000000000000000..d693b74a1945e17f4ef454030551804b5f953b1a
--- /dev/null
+++ b/regtest/basic/rt-mpi7/rgrid.0
@@ -0,0 +1,38 @@
+#! FIELDS d1 d2 restart.bias der_d1 der_d2
+#! SET min_d1 0
+#! SET max_d1 2
+#! SET nbins_d1 5
+#! SET periodic_d1 false
+#! SET min_d2 0
+#! SET max_d2 2
+#! SET nbins_d2 5
+#! SET periodic_d2 false
+ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.000000000 0.002715284 0.086233767 0.015955653
+ 1.500000000 1.000000000 0.435779639 -3.722399196 5.324145962
+ 2.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 1.500000000 0.010483058 -0.012211947 -0.306461908
+ 2.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
diff --git a/regtest/basic/rt-mpi7/rgrid.1 b/regtest/basic/rt-mpi7/rgrid.1
new file mode 100644
index 0000000000000000000000000000000000000000..a8975ad3c11dd71847f5f00f6e7832c02538e4a2
--- /dev/null
+++ b/regtest/basic/rt-mpi7/rgrid.1
@@ -0,0 +1,38 @@
+#! FIELDS d1 d2 restart.bias der_d1 der_d2
+#! SET min_d1 0
+#! SET max_d1 2
+#! SET nbins_d1 5
+#! SET periodic_d1 false
+#! SET min_d2 0
+#! SET max_d2 2
+#! SET nbins_d2 5
+#! SET periodic_d2 false
+ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.000000000 0.012633933 0.346690534 0.112369568
+ 1.500000000 1.000000000 0.406125472 -4.119753516 4.462299304
+ 2.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 1.500000000 0.001650351 -0.004046884 -0.055641830
+ 2.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
diff --git a/regtest/basic/rt-mpi7/rgrid.2 b/regtest/basic/rt-mpi7/rgrid.2
new file mode 100644
index 0000000000000000000000000000000000000000..9decb0c312aeb6b7660d032493493ddead04ee42
--- /dev/null
+++ b/regtest/basic/rt-mpi7/rgrid.2
@@ -0,0 +1,38 @@
+#! FIELDS d1 d2 restart.bias der_d1 der_d2
+#! SET min_d1 0
+#! SET max_d1 2
+#! SET nbins_d1 5
+#! SET periodic_d1 false
+#! SET min_d2 0
+#! SET max_d2 2
+#! SET nbins_d2 5
+#! SET periodic_d2 false
+ 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 0.500000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 0.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.000000000 0.009918649 0.260456766 0.096413915
+ 1.500000000 1.000000000 0.374700112 -2.710607032 5.036795586
+ 2.000000000 1.000000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 1.500000000 0.010483058 -0.012211947 -0.306461908
+ 2.000000000 1.500000000 0.000000000 0.000000000 0.000000000
+
+ 0.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 0.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000
+ 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000
diff --git a/scripts/vim2html.sh b/scripts/vim2html.sh
index e8e891bf2a0d0f9cf8010be40db67135a94008c6..d39c2c23aec3a8e5d909e4cf2c474158e26dd290 100755
--- a/scripts/vim2html.sh
+++ b/scripts/vim2html.sh
@@ -18,7 +18,10 @@ Usage:
Options can be:
--annotate-syntax Reports annotated syntax on output.
Also reports non-annotated regions on stderr
-
+--pdf To produce a pdf file with syntax highlighting.
+--crop Crop the pdf file to the only written part. Usefull to insert the pdf in a LaTex file as image.
+--fs Specify the fontsize of the pdf output.
+--colors Specify the color palette. Allowed values are: default/ac
EOF
exit 0
fi
@@ -31,12 +34,30 @@ fi
annotate_syntax=no
inputfile=""
outputfile=""
-
+pdf=no
+crop=no
+fontsize=7.5
+prefix=""
+colors=""
for opt
do
- case "$opt" in
+ prefixopt="$prefix$opt"
+ prefix=""
+ case "$prefixopt" in
(--annotate-syntax)
annotate_syntax=yes ;;
+ (--pdf)
+ pdf=yes;;
+ (--crop)
+ crop=yes;;
+ (--fs)
+ prefix="--fs=";;
+ (--fs=*)
+ fontsize="${prefixopt#--fs=}";;
+ (--colors)
+ prefix="--colors=";;
+ (--colors=*)
+ colors="${prefixopt#--colors=}";;
(-*)
echo "ERROR: Unknown option $opt. Use -h for help."
exit 1 ;;
@@ -51,7 +72,6 @@ do
fi
esac
done
-
temp="$( mktemp -dt plumed.XXXXX)"
{
@@ -64,45 +84,95 @@ fi | tr '\r' ' '
} | {
cd $temp
{
-if [ "$annotate_syntax" = "yes" ]; then
-vim - -c ":syntax off" \
+case "$colors" in
+(ac)
+ cat > $temp/colors.vim << EOF
+: set bg=dark
+: hi Statement ctermfg=5 cterm=bold guifg=#ffff60 gui=none
+: highlight Type ctermfg=2 cterm=bold
+: hi Comment guifg=#80a0ff ctermfg=darkblue
+: hi Constant ctermfg=red guifg=#ffa0a0 cterm=bold
+EOF
+;;
+(*)
+ echo > $temp/colors.vim
+;;
+esac
+
+if [ "$pdf" = "yes" ]; then
+
+# this is a workaround for a bug in TOhtml
+# that does not color properly when there are empty lines
+ sed "s/^$/ /" |
+ vim - -c ":syntax off" \
+ -c ":source $VIMFILE" \
+ -c ":source $temp/colors.vim" \
+ -c ":set printoptions=header:0" \
+ -c ":set printfont=courier:h$fontsize" \
+ -c ":hardcopy > out.ps" \
+ -c ":q!" 1> /dev/null 2> /dev/null
+else
+ if [ "$annotate_syntax" = "yes" ]; then
+
+ vim - -c ":syntax off" \
-c ":source $VIMFILE" \
-c ":call PlumedAnnotateSyntax()" \
-c ":wq! syntax.log" -c ":q!" 1> /dev/null 2> /dev/null
-else
-# this is a workaround for a bug in TOhtml
-# that does not color properly when there are empty lines
-sed "s/^$/ /" |
-vim - -c ":syntax off" \
+ else
+ # this is a workaround for a bug in TOhtml
+ # that does not color properly when there are empty lines
+ sed "s/^$/ /" |
+ vim - -c ":syntax off" \
-c ":source $VIMFILE" \
+ -c ":source $temp/colors.vim" \
-c ":let g:html_use_css = 0" \
-c ":let g:html_no_progress" \
-c ":let g:html_number_lines = 1" \
-c ":TOhtml" \
-c ":wq! out.html" -c ":q!" 1> /dev/null 2> /dev/null
+ fi
fi
+
+
}
}
-if [ "$annotate_syntax" = "yes" ] ; then
+if [ "$pdf" = "yes" ]; then
+ ps2pdf $temp/out.ps $temp/out.pdf
+ output="$temp/out.pdf"
+ if [ "$crop" = "yes" ]; then
+ pdfcrop $output $temp/outcrop.pdf
+ output="$temp/outcrop.pdf"
+ fi
+else
+ if [ "$annotate_syntax" = "yes" ] ; then
+ output="$(awk '{if($1=="ANNOTATION") print}' $temp/syntax.log 2>/dev/null)"
+ error="$(awk '{if($1=="ERROR") print}' $temp/syntax.log >&2 2>/dev/null)"
+ else
+ output="$(cat $temp/out.html | sed 's|^<title>.*</title>$||' | sed 's|^<meta.*$||')"
+ fi
+fi
-output="$(awk '{if($1=="ANNOTATION") print}' $temp/syntax.log 2>/dev/null)"
-error="$(awk '{if($1=="ERROR") print}' $temp/syntax.log >&2 2>/dev/null)"
+if [ "$pdf" = "yes" ]; then
+ if [ -n "$outputfile" ] ; then
+ cp $output $outputfile
+ else
+ echo "You should specify an output file if using pdf opton!"
+ exit 1;
+ fi
else
-output="$(cat $temp/out.html | sed 's|^<title>.*</title>$||' | sed 's|^<meta.*$||')"
+ if [ -n "$outputfile" ] ; then
+ echo "$output" > "$outputfile"
+ else
+ echo "$output"
+ fi
fi
-if [ -n "$outputfile" ] ; then
- echo "$output" > "$outputfile"
-else
- echo "$output"
-fi
-
echo "$error" >&2
if [ -n "$error" ] ; then
diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index d6c2c1a6c52608e82b804b60a0aa278e1279e6e4..eae4272b81459803a45cb8f9f24de663e39eaf0e 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -85,6 +85,10 @@ In case you do not provide any information about bin size (neither GRID_BIN nor
and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing.
This default choice should be reasonable for most applications.
+Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second
+case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read
+it using GRID_RFILE.
+
Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this
varient of metadynamics the heights of the Gaussian hills are rescaled at each step so the bias is now
given by:
@@ -272,7 +276,7 @@ private:
Gaussian(const vector<double> & center,const vector<double> & sigma,double height, bool multivariate ):
center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma){
// to avoid troubles from zero element in flexible hills
- for(unsigned i=0;i<invsigma.size();++i)abs(invsigma[i])>1.e-20?invsigma[i]=1.0/invsigma[i]:0.;
+ for(unsigned i=0;i<invsigma.size();++i) abs(invsigma[i])>1.e-20?invsigma[i]=1.0/invsigma[i]:0.;
}
};
vector<double> sigma0_;
@@ -282,11 +286,9 @@ private:
OFile hillsOfile_;
OFile gridfile_;
Grid* BiasGrid_;
- Grid* ExtGrid_;
bool storeOldGrids_;
- std::string gridfilename_,gridreadfilename_;
int wgridstride_;
- bool grid_,hasextgrid_;
+ bool grid_;
double height0_;
double biasf_;
double kbt_;
@@ -309,12 +311,10 @@ private:
bool doInt_;
bool isFirstStep;
double reweight_factor;
- std::vector<unsigned> rewf_grid_;
- unsigned int rewf_ustride_;
-/// accumulator for work
+ vector<unsigned> rewf_grid_;
+ unsigned rewf_ustride_;
double work_;
long int last_step_warn_grid;
-
void readGaussians(IFile*);
bool readChunkOfGaussians(IFile *ifile, unsigned n);
@@ -326,9 +326,9 @@ private:
void finiteDifferenceGaussian(const vector<double>&, const Gaussian&);
double getGaussianNormalization( const Gaussian& );
vector<unsigned> getGaussianSupport(const Gaussian&);
- bool scanOneHill(IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate );
+ bool scanOneHill(IFile *ifile, vector<Value> &v, vector<double> ¢er, vector<double> &sigma, double &height, bool &multivariate);
void computeReweightingFactor();
- std::string fmt;
+ string fmt;
public:
explicit MetaD(const ActionOptions&);
@@ -345,7 +345,8 @@ void MetaD::registerKeywords(Keywords& keys){
Bias::registerKeywords(keys);
componentsAreNotOptional(keys);
keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
- keys.addOutputComponent("rbias","REWEIGHTING_NGRID","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-c(t)]. This component can be used to obtain a reweighted histogram.");
+ keys.addOutputComponent("rbias","REWEIGHTING_NGRID","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-c(t)]."
+ "This component can be used to obtain a reweighted histogram.");
keys.addOutputComponent("rct","REWEIGHTING_NGRID","the reweighting factor \\f$c(t)\\f$.");
keys.addOutputComponent("work","default","accumulator for work");
keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor");
@@ -362,12 +363,17 @@ void MetaD::registerKeywords(Keywords& keys){
keys.add("optional","GRID_MAX","the upper bounds for the grid");
keys.add("optional","GRID_BIN","the number of bins for the grid");
keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
- keys.add("optional","REWEIGHTING_NGRID","calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-c(t)]. Here you should specify the number of grid points required in each dimension. The number of grid points should be equal or larger to the number of grid points given in GRID_BIN." "This method is not compatible with metadynamics not on a grid.");
- keys.add("optional","REWEIGHTING_NHILLS","how many Gaussian hills should be deposited between calculating the c(t) reweighting factor. The default is to do this every 50 hills.");
+ keys.add("optional","REWEIGHTING_NGRID","calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-c(t)]."
+ "Here you should specify the number of grid points required in each dimension."
+ "The number of grid points should be equal or larger to the number of grid points given in GRID_BIN."
+ "This method is not compatible with metadynamics not on a grid.");
+ keys.add("optional","REWEIGHTING_NHILLS","how many Gaussian hills should be deposited between calculating the c(t) reweighting factor."
+ "The default is to do this every 50 hills.");
keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills");
keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids");
keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps");
keys.add("optional","GRID_WFILE","the file on which to write the grid");
+ keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present");
keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or timestep dimensions");
keys.add("optional","WALKERS_ID", "walker id");
@@ -375,7 +381,6 @@ void MetaD::registerKeywords(Keywords& keys){
keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers");
keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files");
keys.add("optional","INTERVAL","monodimensional lower and upper limits, outside the limits the system will not feel the biasing force.");
- keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation");
keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds ");
keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with other WALKERS_* options");
@@ -401,7 +406,7 @@ MetaD::~MetaD(){
MetaD::MetaD(const ActionOptions& ao):
PLUMED_BIAS_INIT(ao),
// Grid stuff initialization
-BiasGrid_(NULL),ExtGrid_(NULL), wgridstride_(0), grid_(false), hasextgrid_(false),
+BiasGrid_(NULL), wgridstride_(0), grid_(false),
// Metadynamics basic parameters
height0_(std::numeric_limits<double>::max()), biasf_(1.0), kbt_(0.0),
stride_(0), welltemp_(false),
@@ -424,55 +429,54 @@ last_step_warn_grid(0)
string adaptiveoption;
adaptiveoption="NONE";
parse("ADAPTIVE",adaptiveoption);
- if (adaptiveoption=="GEOM"){
- log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
- adaptive_=FlexibleBin::geometry;
- }else if (adaptiveoption=="DIFF"){
- log.printf(" Uses Diffusion-based hills width: sigma must be in timesteps and only one sigma is needed\n");
- adaptive_=FlexibleBin::diffusion;
- }else if (adaptiveoption=="NONE"){
- adaptive_=FlexibleBin::none;
- }else{
- error("I do not know this type of adaptive scheme");
+ if(adaptiveoption=="GEOM"){
+ log.printf(" Uses Geometry-based hills width: sigma must be in distance units and only one sigma is needed\n");
+ adaptive_=FlexibleBin::geometry;
+ } else if(adaptiveoption=="DIFF"){
+ log.printf(" Uses Diffusion-based hills width: sigma must be in timesteps and only one sigma is needed\n");
+ adaptive_=FlexibleBin::diffusion;
+ } else if(adaptiveoption=="NONE"){
+ adaptive_=FlexibleBin::none;
+ } else {
+ error("I do not know this type of adaptive scheme");
}
// parse the sigma
parseVector("SIGMA",sigma0_);
parse("FMT",fmt);
- if (adaptive_==FlexibleBin::none){
- // if you use normal sigma you need one sigma per argument
- if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
- }else{
- // if you use flexible hills you need one sigma
- if(sigma0_.size()!=1){
- error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
- }
- // if adaptive then the number must be an integer
- if(adaptive_==FlexibleBin::diffusion){
- if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ){
- error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of timesteps\n");
- }
- }
- // here evtl parse the sigma min and max values
-
- parseVector("SIGMA_MIN",sigma0min_);
- if(sigma0min_.size()>0 && sigma0min_.size()<getNumberOfArguments()) {
- error("the number of SIGMA_MIN values be at least the number of the arguments");
- } else if(sigma0min_.size()==0) {
- sigma0min_.resize(getNumberOfArguments());
- for(unsigned i=0;i<getNumberOfArguments();i++){sigma0min_[i]=-1.;}
- }
-
- parseVector("SIGMA_MAX",sigma0max_);
- if(sigma0max_.size()>0 && sigma0max_.size()<getNumberOfArguments()) {
- error("the number of SIGMA_MAX values be at least the number of the arguments");
- } else if(sigma0max_.size()==0) {
- sigma0max_.resize(getNumberOfArguments());
- for(unsigned i=0;i<getNumberOfArguments();i++){sigma0max_[i]=-1.;}
- }
-
- flexbin=new FlexibleBin(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_);
+ if(adaptive_==FlexibleBin::none){
+ // if you use normal sigma you need one sigma per argument
+ if( sigma0_.size()!=getNumberOfArguments() ) error("number of arguments does not match number of SIGMA parameters");
+ } else {
+ // if you use flexible hills you need one sigma
+ if(sigma0_.size()!=1){
+ error("If you choose ADAPTIVE you need only one sigma according to your choice of type (GEOM/DIFF)");
+ }
+ // if adaptive then the number must be an integer
+ if(adaptive_==FlexibleBin::diffusion){
+ if(int(sigma0_[0])-sigma0_[0]>1.e-9 || int(sigma0_[0])-sigma0_[0] <-1.e-9 || int(sigma0_[0])<1 ){
+ error("In case of adaptive hills with diffusion, the sigma must be an integer which is the number of timesteps\n");
+ }
+ }
+ // here evtl parse the sigma min and max values
+ parseVector("SIGMA_MIN",sigma0min_);
+ if(sigma0min_.size()>0 && sigma0min_.size()<getNumberOfArguments()) {
+ error("the number of SIGMA_MIN values be at least the number of the arguments");
+ } else if(sigma0min_.size()==0) {
+ sigma0min_.resize(getNumberOfArguments());
+ for(unsigned i=0;i<getNumberOfArguments();i++){sigma0min_[i]=-1.;}
+ }
+
+ parseVector("SIGMA_MAX",sigma0max_);
+ if(sigma0max_.size()>0 && sigma0max_.size()<getNumberOfArguments()) {
+ error("the number of SIGMA_MAX values be at least the number of the arguments");
+ } else if(sigma0max_.size()==0) {
+ sigma0max_.resize(getNumberOfArguments());
+ for(unsigned i=0;i<getNumberOfArguments();i++){sigma0max_[i]=-1.;}
+ }
+
+ flexbin=new FlexibleBin(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_);
}
// note: HEIGHT is not compulsory, since one could use the TAU keyword, see below
parse("HEIGHT",height0_);
@@ -556,6 +560,7 @@ last_step_warn_grid(0)
bool spline=!nospline;
if(gbin.size()>0){grid_=true;}
parse("GRID_WSTRIDE",wgridstride_);
+ string gridfilename_;
parse("GRID_WFILE",gridfilename_);
parseFlag("STORE_GRIDS",storeOldGrids_);
if(grid_ && gridfilename_.length()>0){
@@ -565,20 +570,23 @@ last_step_warn_grid(0)
if(grid_ && wgridstride_>0){
if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE");
}
-
+ string gridreadfilename_;
parse("GRID_RFILE",gridreadfilename_);
+ if(!grid_&&gridfilename_.length()> 0) error("To write a grid you need first to define it!");
+ if(!grid_&&gridreadfilename_.length()>0) error("To read a grid you need first to define it!");
+
if(grid_){
- parseVector("REWEIGHTING_NGRID",rewf_grid_);
- if( rewf_grid_.size()>0 && rewf_grid_.size()!=getNumberOfArguments() ){
- error("size mismatch for REWEIGHTING_NGRID keyword");
- } else if( rewf_grid_.size()==getNumberOfArguments() ){
- for(unsigned j=0;j<getNumberOfArguments();++j){
- if( !getPntrToArgument(j)->isPeriodic() ) rewf_grid_[j] += 1;
- }
- }
- if( adaptive_==FlexibleBin::diffusion || adaptive_==FlexibleBin::geometry) warning("reweighting has not been proven to work with adaptive Gaussians");
- rewf_ustride_=50; parse("REWEIGHTING_NHILLS",rewf_ustride_);
+ parseVector("REWEIGHTING_NGRID",rewf_grid_);
+ if(rewf_grid_.size()>0 && rewf_grid_.size()!=getNumberOfArguments()){
+ error("size mismatch for REWEIGHTING_NGRID keyword");
+ } else if(rewf_grid_.size()==getNumberOfArguments()){
+ for(unsigned j=0;j<getNumberOfArguments();++j){
+ if( !getPntrToArgument(j)->isPeriodic() ) rewf_grid_[j] += 1;
+ }
+ }
+ if(adaptive_==FlexibleBin::diffusion || adaptive_==FlexibleBin::geometry) warning("reweighting has not been proven to work with adaptive Gaussians");
+ rewf_ustride_=50; parse("REWEIGHTING_NHILLS",rewf_ustride_);
}
// Multiple walkers
@@ -622,38 +630,36 @@ last_step_warn_grid(0)
}
if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_);
if(grid_){
- log.printf(" Grid min");
- for(unsigned i=0;i<gmin.size();++i) log.printf(" %s",gmin[i].c_str() );
- log.printf("\n");
- log.printf(" Grid max");
- for(unsigned i=0;i<gmax.size();++i) log.printf(" %s",gmax[i].c_str() );
- log.printf("\n");
- log.printf(" Grid bin");
- for(unsigned i=0;i<gbin.size();++i) log.printf(" %u",gbin[i]);
- log.printf("\n");
- if(spline){log.printf(" Grid uses spline interpolation\n");}
- if(sparsegrid){log.printf(" Grid uses sparse grid\n");}
- if(wgridstride_>0){log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
+ log.printf(" Grid min");
+ for(unsigned i=0;i<gmin.size();++i) log.printf(" %s",gmin[i].c_str() );
+ log.printf("\n");
+ log.printf(" Grid max");
+ for(unsigned i=0;i<gmax.size();++i) log.printf(" %s",gmax[i].c_str() );
+ log.printf("\n");
+ log.printf(" Grid bin");
+ for(unsigned i=0;i<gbin.size();++i) log.printf(" %u",gbin[i]);
+ log.printf("\n");
+ if(spline){log.printf(" Grid uses spline interpolation\n");}
+ if(sparsegrid){log.printf(" Grid uses sparse grid\n");}
+ if(wgridstride_>0){log.printf(" Grid is written on file %s with stride %d\n",gridfilename_.c_str(),wgridstride_);}
}
- if(gridreadfilename_.length()>0){
- log.printf(" Reading an additional bias from grid in file %s \n",gridreadfilename_.c_str());
- }
-
if(mw_n_>1){
- if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
- log.printf(" %d multiple walkers active\n",mw_n_);
- log.printf(" walker id %d\n",mw_id_);
- log.printf(" reading stride %d\n",mw_rstride_);
- log.printf(" directory with hills files %s\n",mw_dir_.c_str());
+ if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers");
+ log.printf(" %d multiple walkers active\n",mw_n_);
+ log.printf(" walker id %d\n",mw_id_);
+ log.printf(" reading stride %d\n",mw_rstride_);
+ log.printf(" directory with hills files %s\n",mw_dir_.c_str());
} else {
- if(walkers_mpi) log.printf(" Multiple walkers active using MPI communnication\n");
+ if(walkers_mpi) log.printf(" Multiple walkers active using MPI communnication\n");
}
addComponent("bias"); componentIsNotPeriodic("bias");
if( rewf_grid_.size()>0 ){
- addComponent("rbias"); componentIsNotPeriodic("rbias");
- addComponent("rct"); componentIsNotPeriodic("rct");
+ addComponent("rbias"); componentIsNotPeriodic("rbias");
+ addComponent("rct"); componentIsNotPeriodic("rct");
+ log.printf(" the c(t) reweighting factor will be calculated every %u hills\n",rewf_ustride_);
+ getPntrToComponent("rct")->set(reweight_factor);
}
addComponent("work"); componentIsNotPeriodic("work");
@@ -663,98 +669,112 @@ last_step_warn_grid(0)
addComponent("acc"); componentIsNotPeriodic("acc");
}
-// for performance
+ // for performance
dp_ = new double[getNumberOfArguments()];
-// initializing and checking grid
- if(grid_){
- // check for adaptive and sigma_min
- if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
- // check for mesh and sigma size
- for(unsigned i=0;i<getNumberOfArguments();i++) {
- double a,b;
- Tools::convert(gmin[i],a);
- Tools::convert(gmax[i],b);
- double mesh=(b-a)/((double)gbin[i]);
- if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
- }
- std::string funcl=getLabel() + ".bias";
- if(!sparsegrid){BiasGrid_=new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
- else{BiasGrid_=new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
- std::vector<std::string> actualmin=BiasGrid_->getMin();
- std::vector<std::string> actualmax=BiasGrid_->getMax();
- for(unsigned i=0;i<getNumberOfArguments();i++){
- if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
- if(gmax[i]!=actualmax[i]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[i]<<" to fit periodicity\n";
- }
- }
-
- if(wgridstride_>0){
- gridfile_.link(*this);
- gridfile_.open(gridfilename_);
- }
-
-// initializing external grid
+ // restart from external grid
+ bool restartedFromGrid=false;
if(gridreadfilename_.length()>0){
- hasextgrid_=true;
- // read the grid in input, find the keys
- IFile gridfile; gridfile.open(gridreadfilename_);
- std::string funcl=getLabel() + ".bias";
- ExtGrid_=Grid::create(funcl,getArguments(),gridfile,false,false,true);
- gridfile.close();
- if(ExtGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
- for(unsigned i=0;i<getNumberOfArguments();++i){
- if( getPntrToArgument(i)->isPeriodic()!=ExtGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
- }
+ // read the grid in input, find the keys
+ IFile gridfile;
+ gridfile.link(*this);
+ if(gridfile.FileExist(gridreadfilename_)){
+ gridfile.open(gridreadfilename_);
+ } else {
+ error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!");
+ }
+ std::string funcl=getLabel() + ".bias";
+ BiasGrid_=Grid::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true);
+ gridfile.close();
+ if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
+ for(unsigned i=0;i<getNumberOfArguments();++i){
+ if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
+ double a, b;
+ Tools::convert(gmin[i],a);
+ Tools::convert(gmax[i],b);
+ double mesh=(b-a)/((double)gbin[i]);
+ if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
+ }
+ log.printf(" Restarting from %s:",gridreadfilename_.c_str());
+ if(getRestart()) restartedFromGrid=true;
}
-// Tiwary-Parrinello reweighting factor
- if(rewf_grid_.size()>0){
- log.printf(" the c(t) reweighting factor will be calculated every %u hills\n",rewf_ustride_);
- getPntrToComponent("rct")->set(reweight_factor);
+ // initializing and checking grid
+ if(grid_&&!(gridreadfilename_.length()>0)){
+ // check for adaptive and sigma_min
+ if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified");
+ // check for mesh and sigma size
+ for(unsigned i=0;i<getNumberOfArguments();i++) {
+ double a,b;
+ Tools::convert(gmin[i],a);
+ Tools::convert(gmax[i],b);
+ double mesh=(b-a)/((double)gbin[i]);
+ if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n";
+ }
+ std::string funcl=getLabel() + ".bias";
+ if(!sparsegrid){BiasGrid_=new Grid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
+ else{BiasGrid_=new SparseGrid(funcl,getArguments(),gmin,gmax,gbin,spline,true);}
+ std::vector<std::string> actualmin=BiasGrid_->getMin();
+ std::vector<std::string> actualmax=BiasGrid_->getMax();
+ for(unsigned i=0;i<getNumberOfArguments();i++){
+ if(gmin[i]!=actualmin[i]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[i]<<" to fit periodicity\n";
+ if(gmax[i]!=actualmax[i]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[i]<<" to fit periodicity\n";
+ }
}
-// creating vector of ifile* for hills reading
-// open all files at the beginning and read Gaussians if restarting
+ // creating vector of ifile* for hills reading
+ // open all files at the beginning and read Gaussians if restarting
for(int i=0;i<mw_n_;++i){
- string fname;
- if(mw_n_>1) {
- stringstream out; out << i;
- fname = mw_dir_+"/"+hillsfname+"."+out.str();
- } else {
- fname = hillsfname;
- }
- IFile *ifile = new IFile();
- ifile->link(*this);
- ifiles.push_back(ifile);
- ifilesnames.push_back(fname);
- if(ifile->FileExist(fname)){
- ifile->open(fname);
- if(getRestart()){
- log.printf(" Restarting from %s:",ifilesnames[i].c_str());
- readGaussians(ifiles[i]);
+ string fname;
+ if(mw_n_>1) {
+ stringstream out; out << i;
+ fname = mw_dir_+"/"+hillsfname+"."+out.str();
+ } else {
+ fname = hillsfname;
+ }
+ IFile *ifile = new IFile();
+ ifile->link(*this);
+ ifiles.push_back(ifile);
+ ifilesnames.push_back(fname);
+ if(ifile->FileExist(fname)){
+ ifile->open(fname);
+ if(getRestart()&&!restartedFromGrid){
+ log.printf(" Restarting from %s:",ifilesnames[i].c_str());
+ readGaussians(ifiles[i]);
+ }
+ ifiles[i]->reset(false);
+ // close only the walker own hills file for later writing
+ if(i==mw_id_) ifiles[i]->close();
}
- ifiles[i]->reset(false);
- // close only the walker own hills file for later writing
- if(i==mw_id_) ifiles[i]->close();
- }
}
comm.Barrier();
-
-// this barrier is needed when using walkers_mpi
-// to be sure that all files have been read before
-// backing them up
-// it should not be used when walkers_mpi is false otherwise
-// it would introduce troubles when using replicas without METAD
-// (e.g. in bias exchange with a neutral replica)
-// see issue #168 on github
+ // this barrier is needed when using walkers_mpi
+ // to be sure that all files have been read before
+ // backing them up
+ // it should not be used when walkers_mpi is false otherwise
+ // it would introduce troubles when using replicas without METAD
+ // (e.g. in bias exchange with a neutral replica)
+ // see issue #168 on github
if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier();
-// Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
- if(plumed.getRestart() && rewf_grid_.size()>0 ){computeReweightingFactor();}
+ // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills
+ if(getRestart() && rewf_grid_.size()>0 ) computeReweightingFactor();
-// open hills file for writing
+ // open grid file for writing
+ if(wgridstride_>0){
+ gridfile_.link(*this);
+ if(walkers_mpi){
+ int r=0;
+ if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
+ comm.Bcast(r,0);
+ if(r>0) gridfilename_="/dev/null";
+ gridfile_.enforceSuffix("");
+ }
+ gridfile_.open(gridfilename_);
+ }
+
+ // open hills file for writing
hillsOfile_.link(*this);
if(walkers_mpi){
int r=0;
@@ -771,11 +791,10 @@ last_step_warn_grid(0)
hillsOfile_.addConstantField("upper_int").printField("upper_int",uppI_);
}
hillsOfile_.setHeavyFlush();
-// output periodicities of variables
+ // output periodicities of variables
for(unsigned i=0;i<getNumberOfArguments();++i) hillsOfile_.setupPrintValue( getPntrToArgument(i) );
bool concurrent=false;
-
const ActionSet&actionSet(plumed.getActionSet());
for(ActionSet::const_iterator p=actionSet.begin();p!=actionSet.end();++p) if(dynamic_cast<MetaD*>(*p)){ concurrent=true; break; }
if(concurrent) log<<" You are using concurrent metadynamics\n";
@@ -801,49 +820,50 @@ last_step_warn_grid(0)
void MetaD::readGaussians(IFile *ifile)
{
- unsigned ncv=getNumberOfArguments();
- vector<double> center(ncv);
- vector<double> sigma(ncv);
- double height;
- int nhills=0;
- bool multivariate=false;
-
- std::vector<Value> tmpvalues;
- for(unsigned j=0;j<getNumberOfArguments();++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
-
- while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
- nhills++;
- if(welltemp_){height*=(biasf_-1.0)/biasf_;}
- addGaussian(Gaussian(center,sigma,height,multivariate));
- }
- log.printf(" %d Gaussians read\n",nhills);
+ unsigned ncv=getNumberOfArguments();
+ vector<double> center(ncv);
+ vector<double> sigma(ncv);
+ double height;
+ int nhills=0;
+ bool multivariate=false;
+
+ std::vector<Value> tmpvalues;
+ for(unsigned j=0;j<getNumberOfArguments();++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
+
+ while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
+ nhills++;
+ if(welltemp_){height*=(biasf_-1.0)/biasf_;}
+ addGaussian(Gaussian(center,sigma,height,multivariate));
+ }
+ log.printf(" %d Gaussians read\n",nhills);
}
bool MetaD::readChunkOfGaussians(IFile *ifile, unsigned n)
{
- unsigned ncv=getNumberOfArguments();
- vector<double> center(ncv);
- vector<double> sigma(ncv);
- double height;
- unsigned nhills=0;
- bool multivariate=false;
- std::vector<Value> tmpvalues;
- for(unsigned j=0;j<getNumberOfArguments();++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
-
- while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
- if(welltemp_){height*=(biasf_-1.0)/biasf_;}
- addGaussian(Gaussian(center,sigma,height,multivariate));
- if(nhills==n){
+ unsigned ncv=getNumberOfArguments();
+ vector<double> center(ncv);
+ vector<double> sigma(ncv);
+ double height;
+ unsigned nhills=0;
+ bool multivariate=false;
+ std::vector<Value> tmpvalues;
+ for(unsigned j=0;j<getNumberOfArguments();++j) tmpvalues.push_back( Value( this, getPntrToArgument(j)->getName(), false ) );
+
+ while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)){;
+ if(welltemp_) height*=(biasf_-1.0)/biasf_;
+ addGaussian(Gaussian(center,sigma,height,multivariate));
+ if(nhills==n){
log.printf(" %u Gaussians read\n",nhills);
return true;
- }
- nhills++;
- }
- log.printf(" %u Gaussians read\n",nhills);
- return false;
+ }
+ nhills++;
+ }
+ log.printf(" %u Gaussians read\n",nhills);
+ return false;
}
-void MetaD::writeGaussian(const Gaussian& hill, OFile&file){
+void MetaD::writeGaussian(const Gaussian& hill, OFile&file)
+{
unsigned ncv=getNumberOfArguments();
file.printField("time",getTimeStep()*getStep());
for(unsigned i=0;i<ncv;++i){
@@ -851,190 +871,183 @@ void MetaD::writeGaussian(const Gaussian& hill, OFile&file){
}
if(hill.multivariate){
hillsOfile_.printField("multivariate","true");
- Matrix<double> mymatrix(ncv,ncv);
- unsigned k=0;
- for(unsigned i=0;i<ncv;i++){
- for(unsigned j=i;j<ncv;j++){
- mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
- k++;
- }
- }
- // invert it
- Matrix<double> invmatrix(ncv,ncv);
- Invert(mymatrix,invmatrix);
- // enforce symmetry
- for(unsigned i=0;i<ncv;i++){
- for(unsigned j=i;j<ncv;j++){
- invmatrix(i,j)=invmatrix(j,i);
- }
- }
-
- // do cholesky so to have a "sigma like" number
- Matrix<double> lower(ncv,ncv);
- cholesky(invmatrix,lower); // now this , in band form , is similar to the sigmas
- // loop in band form
- for (unsigned i=0;i<ncv;i++){
- for (unsigned j=0;j<ncv-i;j++){
- file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
- }
- }
+ Matrix<double> mymatrix(ncv,ncv);
+ unsigned k=0;
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=i;j<ncv;j++){
+ // recompose the full inverse matrix
+ mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
+ k++;
+ }
+ }
+ // invert it
+ Matrix<double> invmatrix(ncv,ncv);
+ Invert(mymatrix,invmatrix);
+ // enforce symmetry
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=i;j<ncv;j++){
+ invmatrix(i,j)=invmatrix(j,i);
+ }
+ }
+
+ // do cholesky so to have a "sigma like" number
+ Matrix<double> lower(ncv,ncv);
+ cholesky(invmatrix,lower);
+ // loop in band form
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=0;j<ncv-i;j++){
+ file.printField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
+ }
+ }
} else {
hillsOfile_.printField("multivariate","false");
for(unsigned i=0;i<ncv;++i)
file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]);
}
double height=hill.height;
- if(welltemp_){height*=biasf_/(biasf_-1.0);}
+ if(welltemp_) height*=biasf_/(biasf_-1.0);
file.printField("height",height).printField("biasf",biasf_);
if(mw_n_>1) file.printField("clock",int(std::time(0)));
file.printField();
}
-void MetaD::addGaussian(const Gaussian& hill){
-
- if(!grid_){hills_.push_back(hill);}
- else{
- unsigned ncv=getNumberOfArguments();
- vector<unsigned> nneighb=getGaussianSupport(hill);
- vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
- vector<double> der(ncv);
- vector<double> xx(ncv);
- if(comm.Get_size()==1){
- for(unsigned i=0;i<neighbors.size();++i){
- Grid::index_t ineigh=neighbors[i];
- for(unsigned j=0;j<ncv;++j){der[j]=0.0;}
- BiasGrid_->getPoint(ineigh,xx);
- double bias=evaluateGaussian(xx,hill,&der[0]);
- BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
- }
- } else {
- unsigned stride=comm.Get_size();
- unsigned rank=comm.Get_rank();
- vector<double> allder(ncv*neighbors.size(),0.0);
- vector<double> allbias(neighbors.size(),0.0);
- for(unsigned i=rank;i<neighbors.size();i+=stride){
- Grid::index_t ineigh=neighbors[i];
- BiasGrid_->getPoint(ineigh,xx);
- allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
- }
- comm.Sum(allbias);
- comm.Sum(allder);
- for(unsigned i=0;i<neighbors.size();++i){
- Grid::index_t ineigh=neighbors[i];
- for(unsigned j=0;j<ncv;++j){der[j]=allder[ncv*i+j];}
- BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
+void MetaD::addGaussian(const Gaussian& hill)
+{
+ if(!grid_) hills_.push_back(hill);
+ else {
+ unsigned ncv=getNumberOfArguments();
+ vector<unsigned> nneighb=getGaussianSupport(hill);
+ vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(hill.center,nneighb);
+ vector<double> der(ncv);
+ vector<double> xx(ncv);
+ if(comm.Get_size()==1){
+ for(unsigned i=0;i<neighbors.size();++i){
+ Grid::index_t ineigh=neighbors[i];
+ for(unsigned j=0;j<ncv;++j) der[j]=0.0;
+ BiasGrid_->getPoint(ineigh,xx);
+ double bias=evaluateGaussian(xx,hill,&der[0]);
+ BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
+ }
+ } else {
+ unsigned stride=comm.Get_size();
+ unsigned rank=comm.Get_rank();
+ vector<double> allder(ncv*neighbors.size(),0.0);
+ vector<double> allbias(neighbors.size(),0.0);
+ for(unsigned i=rank;i<neighbors.size();i+=stride){
+ Grid::index_t ineigh=neighbors[i];
+ BiasGrid_->getPoint(ineigh,xx);
+ allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]);
+ }
+ comm.Sum(allbias);
+ comm.Sum(allder);
+ for(unsigned i=0;i<neighbors.size();++i){
+ Grid::index_t ineigh=neighbors[i];
+ for(unsigned j=0;j<ncv;++j){der[j]=allder[ncv*i+j];}
+ BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
+ }
}
}
- }
}
vector<unsigned> MetaD::getGaussianSupport(const Gaussian& hill)
{
- vector<unsigned> nneigh;
- if(doInt_){
- double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[0];
- if(hill.center[0]+cutoff > uppI_ || hill.center[0]-cutoff < lowI_) {
- // in this case, we updated the entire grid to avoid problems
- return BiasGrid_->getNbin();
- } else {
- nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[0])) );
- return nneigh;
- }
- }
-
- // traditional or flexible hill?
- if(hill.multivariate){
- unsigned ncv=getNumberOfArguments();
- unsigned k=0;
- //log<<"------- GET GAUSSIAN SUPPORT --------\n";
- Matrix<double> mymatrix(ncv,ncv);
- for(unsigned i=0;i<ncv;i++){
- for(unsigned j=i;j<ncv;j++){
- mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
- k++;
- }
- }
- // Reinvert so to have the ellipses
- Matrix<double> myinv(ncv,ncv);
- Invert(mymatrix,myinv);
- Matrix<double> myautovec(ncv,ncv);
- vector<double> myautoval(ncv); //should I take this or their square root?
- diagMat(myinv,myautoval,myautovec);
- double maxautoval;maxautoval=0.;
- unsigned ind_maxautoval;ind_maxautoval=ncv;
- for (unsigned i=0;i<ncv;i++){
- if(myautoval[i]>maxautoval){maxautoval=myautoval[i];ind_maxautoval=i;}
- }
- for (unsigned i=0;i<ncv;i++){
- double cutoff=sqrt(2.0*DP2CUTOFF)*abs(sqrt(maxautoval)*myautovec(i,ind_maxautoval));
- //log<<"AUTOVAL "<<myautoval[0]<<" COMP "<<abs(myautoval[0]*myautovec(i,0)) <<" CUTOFF "<<cutoff<<"\n";
- nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[i])) );
- }
- }else{
- for(unsigned i=0;i<getNumberOfArguments();++i){
- double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[i];
- nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[i])) );
- }
- }
- //log<<"------- END GET GAUSSIAN SUPPORT --------\n";
- return nneigh;
+ vector<unsigned> nneigh;
+ if(doInt_){
+ double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[0];
+ if(hill.center[0]+cutoff > uppI_ || hill.center[0]-cutoff < lowI_) {
+ // in this case, we updated the entire grid to avoid problems
+ return BiasGrid_->getNbin();
+ } else {
+ nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[0])) );
+ return nneigh;
+ }
+ }
+
+ // traditional or flexible hill?
+ if(hill.multivariate){
+ unsigned ncv=getNumberOfArguments();
+ unsigned k=0;
+ Matrix<double> mymatrix(ncv,ncv);
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=i;j<ncv;j++){
+ // recompose the full inverse matrix
+ mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k];
+ k++;
+ }
+ }
+ // Reinvert so to have the ellipses
+ Matrix<double> myinv(ncv,ncv);
+ Invert(mymatrix,myinv);
+ Matrix<double> myautovec(ncv,ncv);
+ vector<double> myautoval(ncv); //should I take this or their square root?
+ diagMat(myinv,myautoval,myautovec);
+ double maxautoval=0.;
+ unsigned ind_maxautoval;ind_maxautoval=ncv;
+ for(unsigned i=0;i<ncv;i++){
+ if(myautoval[i]>maxautoval){maxautoval=myautoval[i];ind_maxautoval=i;}
+ }
+ for(unsigned i=0;i<ncv;i++){
+ const double cutoff=sqrt(2.0*DP2CUTOFF)*abs(sqrt(maxautoval)*myautovec(i,ind_maxautoval));
+ nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[i])) );
+ }
+ } else {
+ for(unsigned i=0;i<getNumberOfArguments();++i){
+ const double cutoff=sqrt(2.0*DP2CUTOFF)*hill.sigma[i];
+ nneigh.push_back( static_cast<unsigned>(ceil(cutoff/BiasGrid_->getDx()[i])) );
+ }
+ }
+ return nneigh;
}
double MetaD::getBiasAndDerivatives(const vector<double>& cv, double* der)
{
- double bias=0.0;
- if(!grid_){
- if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000){
- std::string msg;
- Tools::convert(hills_.size(),msg);
- msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
- warning(msg);
- last_step_warn_grid=getStep();
- }
- unsigned stride=comm.Get_size();
- unsigned rank=comm.Get_rank();
- for(unsigned i=rank;i<hills_.size();i+=stride){
- bias+=evaluateGaussian(cv,hills_[i],der);
- //finite difference test
- //finiteDifferenceGaussian(cv,hills_[i]);
- }
- comm.Sum(bias);
- if(der) comm.Sum(der,getNumberOfArguments());
- }else{
- if(der){
- vector<double> vder(getNumberOfArguments());
- bias=BiasGrid_->getValueAndDerivatives(cv,vder);
- for(unsigned i=0;i<getNumberOfArguments();++i) {der[i]=vder[i];}
- }else{
- bias=BiasGrid_->getValue(cv);
+ double bias=0.0;
+ if(!grid_){
+ if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000){
+ std::string msg;
+ Tools::convert(hills_.size(),msg);
+ msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits";
+ warning(msg);
+ last_step_warn_grid=getStep();
+ }
+ unsigned stride=comm.Get_size();
+ unsigned rank=comm.Get_rank();
+ for(unsigned i=rank;i<hills_.size();i+=stride){
+ bias+=evaluateGaussian(cv,hills_[i],der);
+ //finite difference test
+ //finiteDifferenceGaussian(cv,hills_[i]);
+ }
+ comm.Sum(bias);
+ if(der) comm.Sum(der,getNumberOfArguments());
+ } else {
+ if(der){
+ vector<double> vder(getNumberOfArguments());
+ bias=BiasGrid_->getValueAndDerivatives(cv,vder);
+ for(unsigned i=0;i<getNumberOfArguments();++i) {der[i]=vder[i];}
+ } else {
+ bias = BiasGrid_->getValue(cv);
+ }
}
- }
- if(hasextgrid_){
- if(der){
- vector<double> vder(getNumberOfArguments());
- bias+=ExtGrid_->getValueAndDerivatives(cv,vder);
- for(unsigned i=0;i<getNumberOfArguments();++i) {der[i]+=vder[i];}
- }else{
- bias+=ExtGrid_->getValue(cv);
- }
- }
- return bias;
+
+ return bias;
}
-double MetaD::getGaussianNormalization( const Gaussian& hill ){
- double norm=1; unsigned ncv=hill.center.size();
+double MetaD::getGaussianNormalization( const Gaussian& hill )
+{
+ double norm=1;
+ unsigned ncv=hill.center.size();
if(hill.multivariate){
- // recompose the full sigma from the upper diag cholesky
+ // recompose the full sigma from the upper diag cholesky
unsigned k=0;
Matrix<double> mymatrix(ncv,ncv);
for(unsigned i=0;i<ncv;i++){
- for(unsigned j=i;j<ncv;j++){
- mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
- k++;
- }
- double ldet; logdet( mymatrix, ldet );
- norm = exp( ldet ); // Not sure here if mymatrix is sigma or inverse
+ for(unsigned j=i;j<ncv;j++){
+ mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
+ k++;
+ }
+ double ldet; logdet( mymatrix, ldet );
+ norm = exp( ldet ); // Not sure here if mymatrix is sigma or inverse
}
} else {
for(unsigned i=0;i<hill.sigma.size();++i) norm*=hill.sigma[i];
@@ -1043,120 +1056,119 @@ double MetaD::getGaussianNormalization( const Gaussian& hill ){
return norm*pow(2*pi,static_cast<double>(ncv)/2.0);
}
-double MetaD::evaluateGaussian
- (const vector<double>& cv, const Gaussian& hill, double* der)
+double MetaD::evaluateGaussian(const vector<double>& cv, const Gaussian& hill, double* der)
{
- double dp2=0.0;
- double bias=0.0;
-// I use a pointer here because cv is const (and should be const)
-// but when using doInt it is easier to locally replace cv[0] with
-// the upper/lower limit in case it is out of range
- const double *pcv=NULL; // pointer to cv
- double tmpcv[1]; // tmp array with cv (to be used with doInt_)
- if(cv.size()>0) pcv=&cv[0];
- if(doInt_){
- plumed_assert(cv.size()==1);
- tmpcv[0]=cv[0];
- if(cv[0]<lowI_) tmpcv[0]=lowI_;
- if(cv[0]>uppI_) tmpcv[0]=uppI_;
- pcv=&(tmpcv[0]);
- }
- if(hill.multivariate){
+ double dp2=0.0;
+ double bias=0.0;
+ // I use a pointer here because cv is const (and should be const)
+ // but when using doInt it is easier to locally replace cv[0] with
+ // the upper/lower limit in case it is out of range
+ const double *pcv=NULL; // pointer to cv
+ double tmpcv[1]; // tmp array with cv (to be used with doInt_)
+ if(cv.size()>0) pcv=&cv[0];
+ if(doInt_){
+ plumed_assert(cv.size()==1);
+ tmpcv[0]=cv[0];
+ if(cv[0]<lowI_) tmpcv[0]=lowI_;
+ if(cv[0]>uppI_) tmpcv[0]=uppI_;
+ pcv=&(tmpcv[0]);
+ }
+ if(hill.multivariate){
unsigned k=0;
unsigned ncv=cv.size();
// recompose the full sigma from the upper diag cholesky
Matrix<double> mymatrix(ncv,ncv);
for(unsigned i=0;i<ncv;i++){
- for(unsigned j=i;j<ncv;j++){
- mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
- k++;
- }
+ for(unsigned j=i;j<ncv;j++){
+ mymatrix(i,j)=mymatrix(j,i)=hill.sigma[k]; // recompose the full inverse matrix
+ k++;
+ }
}
-
for(unsigned i=0;i<cv.size();++i){
- double dp_i=difference(i,hill.center[i],pcv[i]);
- dp_[i]=dp_i;
- for(unsigned j=i;j<cv.size();++j){
- if(i==j){
- dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
- }else{
- double dp_j=difference(j,hill.center[j],pcv[j]);
- dp2+=dp_i*dp_j*mymatrix(i,j) ;
- }
+ double dp_i=difference(i,hill.center[i],pcv[i]);
+ dp_[i]=dp_i;
+ for(unsigned j=i;j<cv.size();++j){
+ if(i==j){
+ dp2+=dp_i*dp_i*mymatrix(i,j)*0.5;
+ } else {
+ double dp_j=difference(j,hill.center[j],pcv[j]);
+ dp2+=dp_i*dp_j*mymatrix(i,j);
}
+ }
}
if(dp2<DP2CUTOFF){
- bias=hill.height*exp(-dp2);
- if(der){
+ bias=hill.height*exp(-dp2);
+ if(der){
for(unsigned i=0;i<cv.size();++i){
- double tmp=0.0;
- k=i;
- for(unsigned j=0;j<cv.size();++j){
- tmp+= dp_[j]*mymatrix(i,j)*bias;
- }
- der[i]-=tmp;
- }
- }
+ double tmp=0.0;
+ k=i;
+ for(unsigned j=0;j<cv.size();++j){
+ tmp += dp_[j]*mymatrix(i,j)*bias;
+ }
+ der[i]-=tmp;
+ }
+ }
}
- }else{
+ } else {
for(unsigned i=0;i<cv.size();++i){
- double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
- dp2+=dp*dp;
- dp_[i]=dp;
+ double dp=difference(i,hill.center[i],pcv[i])*hill.invsigma[i];
+ dp2+=dp*dp;
+ dp_[i]=dp;
}
dp2*=0.5;
if(dp2<DP2CUTOFF){
- bias=hill.height*exp(-dp2);
- if(der){
+ bias=hill.height*exp(-dp2);
+ if(der){
for(unsigned i=0;i<cv.size();++i){der[i]+=-bias*dp_[i]*hill.invsigma[i];}
- }
+ }
}
- }
- if(doInt_){
- if((cv[0]<lowI_ || cv[0]>uppI_) && der ) for(unsigned i=0;i<cv.size();++i)der[i]=0;
}
- return bias;
+
+ if(doInt_ && der){
+ if(cv[0]<lowI_ || cv[0]>uppI_) for(unsigned i=0;i<cv.size();++i) der[i]=0;
+ }
+
+ return bias;
}
double MetaD::getHeight(const vector<double>& cv)
{
- double height=height0_;
- if(welltemp_){
- double vbias=getBiasAndDerivatives(cv);
- height=height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
- }
- return height;
+ double height=height0_;
+ if(welltemp_){
+ double vbias = getBiasAndDerivatives(cv);
+ height = height0_*exp(-vbias/(kbt_*(biasf_-1.0)));
+ }
+ return height;
}
void MetaD::calculate()
{
-
-// this is because presently there is no way to properly pass information
-// on adaptive hills (diff) after exchanges:
+ // this is because presently there is no way to properly pass information
+ // on adaptive hills (diff) after exchanges:
if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange");
- unsigned ncv=getNumberOfArguments();
+ const unsigned ncv=getNumberOfArguments();
vector<double> cv(ncv);
- for(unsigned i=0;i<ncv;++i){cv[i]=getArgument(i);}
-
- double* der=new double[ncv];
- for(unsigned i=0;i<ncv;++i){der[i]=0.0;}
- double ene=getBiasAndDerivatives(cv,der);
+ double* der = new double[ncv];
+ for(unsigned i=0;i<ncv;++i) {
+ cv[i]=getArgument(i);
+ der[i]=0.;
+ }
+ const double ene = getBiasAndDerivatives(cv,der);
getPntrToComponent("bias")->set(ene);
if( rewf_grid_.size()>0 ) getPntrToComponent("rbias")->set(ene - reweight_factor);
-// calculate the acceleration factor
+ // calculate the acceleration factor
if(acceleration&&!isFirstStep) {
acc += exp(ene/(kbt_));
- double mean_acc = acc/((double) getStep());
+ const double mean_acc = acc/((double) getStep());
getPntrToComponent("acc")->set(mean_acc);
}
getPntrToComponent("work")->set(work_);
-// set Forces
+ // set Forces
for(unsigned i=0;i<ncv;++i){
- const double f=-der[i];
- setOutputForce(i,f);
+ const double f=-der[i];
+ setOutputForce(i,f);
}
-
delete [] der;
}
@@ -1167,76 +1179,78 @@ void MetaD::update(){
// adding hills criteria (could be more complex though)
bool nowAddAHill;
- if(getStep()%stride_==0 && !isFirstStep ){nowAddAHill=true;}else{nowAddAHill=false;isFirstStep=false;}
+ if(getStep()%stride_==0 && !isFirstStep )nowAddAHill=true;
+ else {
+ nowAddAHill=false;
+ isFirstStep=false;
+ }
- for(unsigned i=0;i<cv.size();++i){cv[i]=getArgument(i);}
+ for(unsigned i=0;i<cv.size();++i) cv[i] = getArgument(i);
double vbias=getBiasAndDerivatives(cv);
// if you use adaptive, call the FlexibleBin
- if (adaptive_!=FlexibleBin::none){
- flexbin->update(nowAddAHill);
- multivariate=true;
- }else{
- multivariate=false;
- };
+ if(adaptive_!=FlexibleBin::none){
+ flexbin->update(nowAddAHill);
+ multivariate=true;
+ } else {
+ multivariate=false;
+ }
- if(nowAddAHill){ // probably this can be substituted with a signal
- // add a Gaussian
- double height=getHeight(cv);
- // use normal sigma or matrix form?
- if (adaptive_!=FlexibleBin::none){
- thissigma=flexbin->getInverseMatrix(); // returns upper diagonal inverse
- //cerr<<"ADDING HILLS "<<endl;
- }else{
- thissigma=sigma0_; // returns normal sigma
- }
-// In case we use walkers_mpi, it is now necessary to communicate with other replicas.
- if(walkers_mpi){
- int nw=0;
- int mw=0;
- if(comm.Get_rank()==0){
-// Only root of group can communicate with other walkers
- nw=multi_sim_comm.Get_size();
- mw=multi_sim_comm.Get_rank();
- }
-// Communicate to the other members of the same group
-// info abount number of walkers and walker index
- comm.Bcast(nw,0);
- comm.Bcast(mw,0);
-// Allocate arrays to store all walkers hills
- std::vector<double> all_cv(nw*cv.size(),0.0);
- std::vector<double> all_sigma(nw*thissigma.size(),0.0);
- std::vector<double> all_height(nw,0.0);
- std::vector<int> all_multivariate(nw,0);
- if(comm.Get_rank()==0){
-// Communicate (only root)
- multi_sim_comm.Allgather(cv,all_cv);
- multi_sim_comm.Allgather(thissigma,all_sigma);
- multi_sim_comm.Allgather(height,all_height);
- multi_sim_comm.Allgather(int(multivariate),all_multivariate);
- }
-// Share info with group members
- comm.Bcast(all_cv,0);
- comm.Bcast(all_sigma,0);
- comm.Bcast(all_height,0);
- comm.Bcast(all_multivariate,0);
- for(int i=0;i<nw;i++){
-// actually add hills one by one
- std::vector<double> cv_now(cv.size());
- std::vector<double> sigma_now(thissigma.size());
- for(unsigned j=0;j<cv.size();j++) cv_now[j]=all_cv[i*cv.size()+j];
- for(unsigned j=0;j<thissigma.size();j++) sigma_now[j]=all_sigma[i*thissigma.size()+j];
- Gaussian newhill=Gaussian(cv_now,sigma_now,all_height[i],all_multivariate[i]);
- addGaussian(newhill);
- writeGaussian(newhill,hillsOfile_);
- }
- } else {
- Gaussian newhill=Gaussian(cv,thissigma,height,multivariate);
- addGaussian(newhill);
-// print on HILLS file
- writeGaussian(newhill,hillsOfile_);
- }
+ if(nowAddAHill){
+ // add a Gaussian
+ double height=getHeight(cv);
+ // returns upper diagonal inverse
+ if(adaptive_!=FlexibleBin::none) thissigma=flexbin->getInverseMatrix();
+ // returns normal sigma
+ else thissigma=sigma0_;
+
+ // In case we use walkers_mpi, it is now necessary to communicate with other replicas.
+ if(walkers_mpi){
+ int nw=0;
+ int mw=0;
+ if(comm.Get_rank()==0){
+ // Only root of group can communicate with other walkers
+ nw=multi_sim_comm.Get_size();
+ mw=multi_sim_comm.Get_rank();
+ }
+ // Communicate to the other members of the same group
+ // info abount number of walkers and walker index
+ comm.Bcast(nw,0);
+ comm.Bcast(mw,0);
+ // Allocate arrays to store all walkers hills
+ std::vector<double> all_cv(nw*cv.size(),0.0);
+ std::vector<double> all_sigma(nw*thissigma.size(),0.0);
+ std::vector<double> all_height(nw,0.0);
+ std::vector<int> all_multivariate(nw,0);
+ if(comm.Get_rank()==0){
+ // Communicate (only root)
+ multi_sim_comm.Allgather(cv,all_cv);
+ multi_sim_comm.Allgather(thissigma,all_sigma);
+ multi_sim_comm.Allgather(height,all_height);
+ multi_sim_comm.Allgather(int(multivariate),all_multivariate);
+ }
+ // Share info with group members
+ comm.Bcast(all_cv,0);
+ comm.Bcast(all_sigma,0);
+ comm.Bcast(all_height,0);
+ comm.Bcast(all_multivariate,0);
+ for(int i=0;i<nw;i++){
+ // actually add hills one by one
+ std::vector<double> cv_now(cv.size());
+ std::vector<double> sigma_now(thissigma.size());
+ for(unsigned j=0;j<cv.size();j++) cv_now[j]=all_cv[i*cv.size()+j];
+ for(unsigned j=0;j<thissigma.size();j++) sigma_now[j]=all_sigma[i*thissigma.size()+j];
+ Gaussian newhill=Gaussian(cv_now,sigma_now,all_height[i],all_multivariate[i]);
+ addGaussian(newhill);
+ writeGaussian(newhill,hillsOfile_);
+ }
+ } else {
+ Gaussian newhill=Gaussian(cv,thissigma,height,multivariate);
+ addGaussian(newhill);
+ // print on HILLS file
+ writeGaussian(newhill,hillsOfile_);
+ }
}
double vbias1=getBiasAndDerivatives(cv);
@@ -1261,167 +1275,159 @@ void MetaD::update(){
}
// if multiple walkers and time to read Gaussians
- if(mw_n_>1 && getStep()%mw_rstride_==0){
- for(int i=0;i<mw_n_;++i){
- // don't read your own Gaussians
- if(i==mw_id_) continue;
- // if the file is not open yet
- if(!(ifiles[i]->isOpen())){
- // check if it exists now and open it!
- if(ifiles[i]->FileExist(ifilesnames[i])) {
- ifiles[i]->open(ifilesnames[i]);
- ifiles[i]->reset(false);
- }
- // otherwise read the new Gaussians
- } else {
- log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
- readGaussians(ifiles[i]);
- ifiles[i]->reset(false);
+ if(mw_n_>1 && getStep()%mw_rstride_==0){
+ for(int i=0;i<mw_n_;++i){
+ // don't read your own Gaussians
+ if(i==mw_id_) continue;
+ // if the file is not open yet
+ if(!(ifiles[i]->isOpen())){
+ // check if it exists now and open it!
+ if(ifiles[i]->FileExist(ifilesnames[i])) {
+ ifiles[i]->open(ifilesnames[i]);
+ ifiles[i]->reset(false);
+ }
+ // otherwise read the new Gaussians
+ } else {
+ log.printf(" Reading hills from %s:",ifilesnames[i].c_str());
+ readGaussians(ifiles[i]);
+ ifiles[i]->reset(false);
+ }
}
- }
- }
-
- if(getStep()%(stride_*rewf_ustride_)==0 && nowAddAHill && rewf_grid_.size()>0 ){computeReweightingFactor();}
+ }
+ if(getStep()%(stride_*rewf_ustride_)==0 && nowAddAHill && rewf_grid_.size()>0 ) computeReweightingFactor();
}
-void MetaD::finiteDifferenceGaussian
- (const vector<double>& cv, const Gaussian& hill)
+void MetaD::finiteDifferenceGaussian(const vector<double>& cv, const Gaussian& hill)
{
- log<<"--------- finiteDifferenceGaussian: size "<<cv.size() <<"------------\n";
- // for each cv
- // first get the bias and the derivative
- vector<double> oldder(cv.size());
- vector<double> der(cv.size());
- vector<double> mycv(cv.size());
- mycv=cv;
- double step=1.e-6;
- Random random;
- // just displace a tiny bit
- for(unsigned i=0;i<cv.size();i++)log<<"CV "<<i<<" V "<<mycv[i]<<"\n";
- for(unsigned i=0;i<cv.size();i++)mycv[i]+=1.e-2*2*(random.RandU01()-0.5);
- for(unsigned i=0;i<cv.size();i++)log<<"NENEWWCV "<<i<<" V "<<mycv[i]<<"\n";
- double oldbias=evaluateGaussian(mycv,hill,&oldder[0]);
- for (unsigned i=0;i<mycv.size();i++){
- double delta=step*2*(random.RandU01()-0.5);
- mycv[i]+=delta;
- double newbias=evaluateGaussian(mycv,hill,&der[0]);
- log<<"CV "<<i;
- log<<" ANAL "<<oldder[i]<<" NUM "<<(newbias-oldbias)/delta<<" DIFF "<<(oldder[i]-(newbias-oldbias)/delta)<<"\n";
- mycv[i]-=delta;
- }
- log<<"--------- END finiteDifferenceGaussian ------------\n";
+ log<<"--------- finiteDifferenceGaussian: size "<<cv.size() <<"------------\n";
+ // for each cv
+ // first get the bias and the derivative
+ vector<double> oldder(cv.size());
+ vector<double> der(cv.size());
+ vector<double> mycv(cv.size());
+ mycv=cv;
+ double step=1.e-6;
+ Random random;
+ // just displace a tiny bit
+ for(unsigned i=0;i<cv.size();i++) log<<"CV "<<i<<" V "<<mycv[i]<<"\n";
+ for(unsigned i=0;i<cv.size();i++) mycv[i]+=1.e-2*2*(random.RandU01()-0.5);
+ for(unsigned i=0;i<cv.size();i++) log<<"NENEWWCV "<<i<<" V "<<mycv[i]<<"\n";
+ double oldbias=evaluateGaussian(mycv,hill,&oldder[0]);
+ for(unsigned i=0;i<mycv.size();i++){
+ double delta=step*2*(random.RandU01()-0.5);
+ mycv[i]+=delta;
+ double newbias=evaluateGaussian(mycv,hill,&der[0]);
+ log<<"CV "<<i;
+ log<<" ANAL "<<oldder[i]<<" NUM "<<(newbias-oldbias)/delta<<" DIFF "<<(oldder[i]-(newbias-oldbias)/delta)<<"\n";
+ mycv[i]-=delta;
+ }
+ log<<"--------- END finiteDifferenceGaussian ------------\n";
}
/// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
-bool MetaD::scanOneHill(IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height , bool &multivariate ){
+bool MetaD::scanOneHill(IFile *ifile, vector<Value> &tmpvalues, vector<double> ¢er, vector<double> &sigma, double &height , bool &multivariate)
+{
double dummy;
multivariate=false;
if(ifile->scanField("time",dummy)){
- unsigned ncv; ncv=tmpvalues.size();
- for(unsigned i=0;i<ncv;++i){
- ifile->scanField( &tmpvalues[i] );
- if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ){
+ unsigned ncv; ncv=tmpvalues.size();
+ for(unsigned i=0;i<ncv;++i){
+ ifile->scanField( &tmpvalues[i] );
+ if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ){
+ error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
+ } else if( tmpvalues[i].isPeriodic() ){
+ std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
+ std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
+ if( imin!=rmin || imax!=rmax ){
error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
- } else if( tmpvalues[i].isPeriodic() ){
- std::string imin, imax; tmpvalues[i].getDomain( imin, imax );
- std::string rmin, rmax; getPntrToArgument(i)->getDomain( rmin, rmax );
- if( imin!=rmin || imax!=rmax ){
- error("in hills file periodicity for variable " + tmpvalues[i].getName() + " does not match periodicity in input");
- }
- }
- center[i]=tmpvalues[i].get();
- }
- // scan for multivariate label: record the actual file position so to eventually rewind
- std::string sss;
- ifile->scanField("multivariate",sss);
- if(sss=="true") multivariate=true;
- else if(sss=="false") multivariate=false;
- else plumed_merror("cannot parse multivariate = "+ sss);
- if(multivariate){
- sigma.resize(ncv*(ncv+1)/2);
- Matrix<double> upper(ncv,ncv);
- Matrix<double> lower(ncv,ncv);
- for (unsigned i=0;i<ncv;i++){
- for (unsigned j=0;j<ncv-i;j++){
- ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
- upper(j,j+i)=lower(j+i,j);
- }
}
- Matrix<double> mymult(ncv,ncv);
- Matrix<double> invmatrix(ncv,ncv);
- //log<<"Lower \n";
- //matrixOut(log,lower);
- //log<<"Upper \n";
- //matrixOut(log,upper);
- mult(lower,upper,mymult);
- //log<<"Mult \n";
- //matrixOut(log,mymult);
- // now invert and get the sigmas
- Invert(mymult,invmatrix);
- //log<<"Invert \n";
- //matrixOut(log,invmatrix);
- // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
- unsigned k=0;
- for (unsigned i=0;i<ncv;i++){
- for (unsigned j=i;j<ncv;j++){
- sigma[k]=invmatrix(i,j);
- k++;
- }
- }
- }else{
- for(unsigned i=0;i<ncv;++i){
- ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
+ }
+ center[i]=tmpvalues[i].get();
+ }
+ // scan for multivariate label: record the actual file position so to eventually rewind
+ std::string sss;
+ ifile->scanField("multivariate",sss);
+ if(sss=="true") multivariate=true;
+ else if(sss=="false") multivariate=false;
+ else plumed_merror("cannot parse multivariate = "+ sss);
+ if(multivariate){
+ sigma.resize(ncv*(ncv+1)/2);
+ Matrix<double> upper(ncv,ncv);
+ Matrix<double> lower(ncv,ncv);
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=0;j<ncv-i;j++){
+ ifile->scanField("sigma_"+getPntrToArgument(j+i)->getName()+"_"+getPntrToArgument(j)->getName(),lower(j+i,j));
+ upper(j,j+i)=lower(j+i,j);
}
- }
-
- ifile->scanField("height",height);
- ifile->scanField("biasf",dummy);
- if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
- if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
- if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
- ifile->scanField();
- return true;
- }else{
+ }
+ Matrix<double> mymult(ncv,ncv);
+ Matrix<double> invmatrix(ncv,ncv);
+ mult(lower,upper,mymult);
+ // now invert and get the sigmas
+ Invert(mymult,invmatrix);
+ // put the sigmas in the usual order: upper diagonal (this time in normal form and not in band form)
+ unsigned k=0;
+ for(unsigned i=0;i<ncv;i++){
+ for(unsigned j=i;j<ncv;j++){
+ sigma[k]=invmatrix(i,j);
+ k++;
+ }
+ }
+ } else {
+ for(unsigned i=0;i<ncv;++i){
+ ifile->scanField("sigma_"+getPntrToArgument(i)->getName(),sigma[i]);
+ }
+ }
+
+ ifile->scanField("height",height);
+ ifile->scanField("biasf",dummy);
+ if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
+ if(ifile->FieldExist("lower_int")) ifile->scanField("lower_int",dummy);
+ if(ifile->FieldExist("upper_int")) ifile->scanField("upper_int",dummy);
+ ifile->scanField();
+ return true;
+ } else {
return false;
- };
+ }
}
-void MetaD::computeReweightingFactor(){
- if( !welltemp_ ) error("cannot compute the c(t) reweighting factors for non well-tempered metadynamics");
+void MetaD::computeReweightingFactor()
+{
+ if( !welltemp_ ) error("cannot compute the c(t) reweighting factors for non well-tempered metadynamics");
- // Recover the minimum values for the grid
- unsigned ncv=getNumberOfArguments();
- unsigned ntotgrid=1;
- std::vector<double> dmin( ncv ),dmax( ncv ), grid_spacing( ncv ), vals( ncv );
- for(unsigned j=0;j<ncv;++j){
+ // Recover the minimum values for the grid
+ unsigned ncv=getNumberOfArguments();
+ unsigned ntotgrid=1;
+ std::vector<double> dmin( ncv ),dmax( ncv ), grid_spacing( ncv ), vals( ncv );
+ for(unsigned j=0;j<ncv;++j){
Tools::convert( BiasGrid_->getMin()[j], dmin[j] );
Tools::convert( BiasGrid_->getMax()[j], dmax[j] );
grid_spacing[j] = ( dmax[j] - dmin[j] ) / static_cast<double>( rewf_grid_[j] );
if( !getPntrToArgument(j)->isPeriodic() ) dmax[j] += grid_spacing[j];
ntotgrid *= rewf_grid_[j];
- }
+ }
- // Now sum over whole grid
- reweight_factor=0.0; double* der=new double[ncv]; std::vector<unsigned> t_index( ncv );
- double sum1=0.0; double sum2=0.0;
- double afactor = biasf_ / (kbt_*(biasf_-1.0)); double afactor2 = 1.0 / (kbt_*(biasf_-1.0));
- unsigned rank=comm.Get_rank(), stride=comm.Get_size();
- for(unsigned i=rank;i<ntotgrid;i+=stride){
- t_index[0]=(i%rewf_grid_[0]);
- unsigned kk=i;
- for(unsigned j=1;j<ncv-1;++j){ kk=(kk-t_index[j-1])/rewf_grid_[i-1]; t_index[j]=(kk%rewf_grid_[i]); }
- if( ncv>=2 ) t_index[ncv-1]=((kk-t_index[ncv-1])/rewf_grid_[ncv-2]);
-
- for(unsigned j=0;j<ncv;++j) vals[j]=dmin[j] + t_index[j]*grid_spacing[j];
-
- double currentb=getBiasAndDerivatives(vals,der);
- sum1 += exp( afactor*currentb );
- sum2 += exp( afactor2*currentb );
- }
- delete [] der;
- comm.Sum( sum1 ); comm.Sum( sum2 );
- reweight_factor = kbt_ * std::log( sum1/sum2 );
- getPntrToComponent("rct")->set(reweight_factor);
+ // Now sum over whole grid
+ reweight_factor=0.0; double* der=new double[ncv]; std::vector<unsigned> t_index( ncv );
+ double sum1=0.0; double sum2=0.0;
+ double afactor = biasf_ / (kbt_*(biasf_-1.0)); double afactor2 = 1.0 / (kbt_*(biasf_-1.0));
+ unsigned rank=comm.Get_rank(), stride=comm.Get_size();
+ for(unsigned i=rank;i<ntotgrid;i+=stride){
+ t_index[0]=(i%rewf_grid_[0]);
+ unsigned kk=i;
+ for(unsigned j=1;j<ncv-1;++j){ kk=(kk-t_index[j-1])/rewf_grid_[i-1]; t_index[j]=(kk%rewf_grid_[i]); }
+ if( ncv>=2 ) t_index[ncv-1]=((kk-t_index[ncv-1])/rewf_grid_[ncv-2]);
+
+ for(unsigned j=0;j<ncv;++j) vals[j]=dmin[j] + t_index[j]*grid_spacing[j];
+
+ double currentb=getBiasAndDerivatives(vals,der);
+ sum1 += exp( afactor*currentb );
+ sum2 += exp( afactor2*currentb );
+ }
+ delete [] der;
+ comm.Sum( sum1 ); comm.Sum( sum2 );
+ reweight_factor = kbt_ * std::log( sum1/sum2 );
+ getPntrToComponent("rct")->set(reweight_factor);
}
}
diff --git a/src/bias/PBMetaD.cpp b/src/bias/PBMetaD.cpp
index b2d045356465d1f26e8d31b0f00070ec35f8f705..b78d9f2e771571c426e08aa9c895bfe7bc353525 100644
--- a/src/bias/PBMetaD.cpp
+++ b/src/bias/PBMetaD.cpp
@@ -177,8 +177,6 @@ MULTIPLE_WALKERS
PRINT ARG=d1,d2,pb.bias STRIDE=100 FILE=COLVAR
\endverbatim
-
->>>>>>> v2.2
*/
//+ENDPLUMEDOC
@@ -195,13 +193,14 @@ private:
vector<double> sigma0_;
vector< vector<Gaussian> > hills_;
vector<OFile*> hillsOfiles_;
- vector<IFile*> ifiles;
+ vector<OFile*> gridfiles_;
vector<Grid*> BiasGrids_;
bool grid_;
double height0_;
double biasf_;
double kbt_;
int stride_;
+ int wgridstride_;
bool welltemp_;
bool multiple_w;
vector<double> uppI_;
@@ -242,6 +241,9 @@ void PBMetaD::registerKeywords(Keywords& keys){
keys.add("optional","BIASFACTOR","use well tempered metadynamics with this biasfactor, one for all biases. Please note you must also specify temp");
keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics");
keys.add("optional","TAU","in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau");
+ keys.add("optional","GRID_RFILES", "read grid for the bias");
+ keys.add("optional","GRID_WFILES", "dump grid for the bias");
+ keys.add("optional","GRID_WSTRIDE", "frequency for dumping the grid");
keys.add("optional","GRID_MIN","the lower bounds for the grid");
keys.add("optional","GRID_MAX","the upper bounds for the grid");
keys.add("optional","GRID_BIN","the number of bins for the grid");
@@ -259,6 +261,12 @@ PBMetaD::~PBMetaD(){
hillsOfiles_[i]->close();
delete hillsOfiles_[i];
}
+ if(wgridstride_ > 0) {
+ for(unsigned i=0; i<gridfiles_.size(); ++i) {
+ gridfiles_[i]->close();
+ delete gridfiles_[i];
+ }
+ }
}
PBMetaD::PBMetaD(const ActionOptions& ao):
@@ -307,6 +315,18 @@ multiple_w(false), doInt_(false), isFirstStep(true)
// MPI version
parseFlag("MULTIPLE_WALKERS",multiple_w);
+
+ // Grid file
+ vector<string> gridfilenames_;
+ parseVector("GRID_WFILES",gridfilenames_);
+ parse("GRID_WSTRIDE",wgridstride_);
+ if (wgridstride_ == 0 && gridfilenames_.size() > 0) {
+ error("frequency with which to output grid not specified use GRID_WSTRIDE");
+ }
+
+ // Read grid
+ vector<string> gridreadfilenames_;
+ parseVector("GRID_RFILES",gridreadfilenames_);
// Grid Stuff
vector<std::string> gmin(getNumberOfArguments());
@@ -326,10 +346,10 @@ multiple_w(false), doInt_(false), isFirstStep(true)
if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present");
if(gmin.size()!=0){
if(gbin.size()==0 && gspacing.size()==0){
- log<<" Binsize not specified, 1/5 of sigma will be be used\n";
+ log<<" Binsize not specified, 1/10 of sigma will be be used\n";
plumed_assert(sigma0_.size()==getNumberOfArguments());
gspacing.resize(getNumberOfArguments());
- for(unsigned i=0;i<gspacing.size();i++) gspacing[i]=0.2*sigma0_[i];
+ for(unsigned i=0;i<gspacing.size();i++) gspacing[i]=0.1*sigma0_[i];
} else if(gspacing.size()!=0 && gbin.size()==0){
log<<" The number of bins will be estimated from GRID_SPACING\n";
} else if(gspacing.size()!=0 && gbin.size()!=0){
@@ -351,7 +371,10 @@ multiple_w(false), doInt_(false), isFirstStep(true)
parseFlag("GRID_NOSPLINE",nospline);
bool spline=!nospline;
if(gbin.size()>0){grid_=true;}
-
+
+ if(!grid_&&gridfilenames_.size() > 0) error("To write a grid you need first to define it!");
+ if(!grid_&&gridreadfilenames_.size() > 0) error("To read a grid you need first to define it!");
+
// Interval keyword
parseVector("INTERVAL_MIN",lowI_);
parseVector("INTERVAL_MAX",uppI_);
@@ -390,14 +413,27 @@ multiple_w(false), doInt_(false), isFirstStep(true)
log.printf("\n");
if(spline){log.printf(" Grid uses spline interpolation\n");}
if(sparsegrid){log.printf(" Grid uses sparse grid\n");}
+ if(wgridstride_>0) {
+ for(unsigned i=0; i<gridfilenames_.size(); ++i) {
+ log.printf(" Grid is written on file %s with stride %d\n",gridfilenames_[i].c_str(),wgridstride_);
+ }
+ }
+ if(gridreadfilenames_.size()>0) {
+ for(unsigned i=0; i<gridreadfilenames_.size(); ++i) {
+ log.printf(" Reading bias from grid in file %s \n",gridreadfilenames_[i].c_str());
+ }
+ }
}
- addComponentWithDerivatives("bias"); componentIsNotPeriodic("bias");
+ addComponent("bias"); componentIsNotPeriodic("bias");
-// initializing vector of hills
+ // initializing vector of hills
hills_.resize(getNumberOfArguments());
-// initializing and checking grid
+ // restart from external grid
+ bool restartedFromGrid=false;
+
+ // initializing and checking grid
if(grid_){
// check for mesh and sigma size
for(unsigned i=0;i<getNumberOfArguments();i++) {
@@ -418,79 +454,122 @@ multiple_w(false), doInt_(false), isFirstStep(true)
gmax_t[0] = gmax[i];
gbin_t[0] = gbin[i];
Grid* BiasGrid_;
- if(!sparsegrid){BiasGrid_=new Grid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
- else {BiasGrid_=new SparseGrid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
- std::vector<std::string> actualmin=BiasGrid_->getMin();
- std::vector<std::string> actualmax=BiasGrid_->getMax();
- if(gmin_t[0]!=actualmin[0]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[0]<<" to fit periodicity\n";
- if(gmax_t[0]!=actualmax[0]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[0]<<" to fit periodicity\n";
+ // Read grid from file
+ if(gridreadfilenames_.size()>0) {
+ IFile gridfile;
+ gridfile.link(*this);
+ if(gridfile.FileExist(gridreadfilenames_[i])){
+ gridfile.open(gridreadfilenames_[i]);
+ } else {
+ error("The GRID file you want to read: " + gridreadfilenames_[i] + ", cannot be found!");
+ }
+ string funcl = getLabel() + ".bias";
+ BiasGrid_ = Grid::create(funcl, args, gridfile, gmin_t, gmax_t, gbin_t, sparsegrid, spline, true);
+ gridfile.close();
+ if(BiasGrid_->getDimension() != args.size()) {
+ error("mismatch between dimensionality of input grid and number of arguments");
+ }
+ if(getPntrToArgument(i)->isPeriodic() != BiasGrid_->getIsPeriodic()[0]) {
+ error("periodicity mismatch between arguments and input bias");
+ }
+ log.printf(" Restarting from %s:",gridreadfilenames_[i].c_str());
+ if(getRestart()) restartedFromGrid=true;
+ } else {
+ if(!sparsegrid){BiasGrid_=new Grid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
+ else {BiasGrid_=new SparseGrid(funcl,args,gmin_t,gmax_t,gbin_t,spline,true);}
+ std::vector<std::string> actualmin=BiasGrid_->getMin();
+ std::vector<std::string> actualmax=BiasGrid_->getMax();
+ if(gmin_t[0]!=actualmin[0]) log<<" WARNING: GRID_MIN["<<i<<"] has been adjusted to "<<actualmin[0]<<" to fit periodicity\n";
+ if(gmax_t[0]!=actualmax[0]) log<<" WARNING: GRID_MAX["<<i<<"] has been adjusted to "<<actualmax[0]<<" to fit periodicity\n";
+ }
BiasGrids_.push_back(BiasGrid_);
}
}
-// read Gaussians if restarting
- for(unsigned i=0;i<hillsfname.size();++i){
- IFile *ifile = new IFile();
- ifile->link(*this);
- if(ifile->FileExist(hillsfname[i])){
- ifile->open(hillsfname[i]);
- if(getRestart()){
- log.printf(" Restarting from %s:",hillsfname[i].c_str());
- readGaussians(i, ifile);
+ // read Gaussians if restarting
+ if (!restartedFromGrid){
+ for(unsigned i=0;i<hillsfname.size();++i){
+ IFile *ifile = new IFile();
+ ifile->link(*this);
+ if(ifile->FileExist(hillsfname[i])){
+ ifile->open(hillsfname[i]);
+ if(getRestart()){
+ log.printf(" Restarting from %s:",hillsfname[i].c_str());
+ readGaussians(i, ifile);
+ }
+ ifile->reset(false);
+ ifile->close();
+ delete ifile;
}
- ifile->reset(false);
- ifile->close();
- delete ifile;
}
}
comm.Barrier();
-
-// this barrier is needed when using walkers_mpi
-// to be sure that all files have been read before
-// backing them up
-// it should not be used when walkers_mpi is false otherwise
-// it would introduce troubles when using replicas without METAD
-// (e.g. in bias exchange with a neutral replica)
-// see issue #168 on github
+ // this barrier is needed when using walkers_mpi
+ // to be sure that all files have been read before
+ // backing them up
+ // it should not be used when walkers_mpi is false otherwise
+ // it would introduce troubles when using replicas without METAD
+ // (e.g. in bias exchange with a neutral replica)
+ // see issue #168 on github
if(comm.Get_rank()==0 && multiple_w) multi_sim_comm.Barrier();
-// open hills files for writing
- for(unsigned i=0;i<hillsfname.size();++i){
- OFile *ofile = new OFile();
- ofile->link(*this);
- string hillsfname_tmp = hillsfname[i];
- // if multiple walkers, only rank 0 will write to file
- if(multiple_w){
+ // open hills files for writing
+ for(unsigned i=0;i<hillsfname.size();++i){
+ OFile *ofile = new OFile();
+ ofile->link(*this);
+ string hillsfname_tmp = hillsfname[i];
+ // if multiple walkers, only rank 0 will write to file
+ if(multiple_w){
int r=0;
if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank();
comm.Bcast(r,0);
if(r>0) hillsfname_tmp="/dev/null";
ofile->enforceSuffix("");
- }
- ofile->open(hillsfname_tmp);
- if(fmt.length()>0) ofile->fmtField(fmt);
- ofile->addConstantField("multivariate");
- if(doInt_) {
+ }
+ ofile->open(hillsfname_tmp);
+ if(fmt.length()>0) ofile->fmtField(fmt);
+ ofile->addConstantField("multivariate");
+ if(doInt_) {
ofile->addConstantField("lower_int").printField("lower_int",lowI_[i]);
ofile->addConstantField("upper_int").printField("upper_int",uppI_[i]);
+ }
+ ofile->setHeavyFlush();
+ // output periodicities of variables
+ ofile->setupPrintValue( getPntrToArgument(i) );
+ // push back
+ hillsOfiles_.push_back(ofile);
+ }
+
+ // Dump grid to files
+ if(wgridstride_ > 0) {
+ for(unsigned i = 0; i < gridfilenames_.size(); ++i) {
+ OFile *ofile = new OFile();
+ ofile->link(*this);
+ string gridfname_tmp = gridfilenames_[i];
+ if(multiple_w) {
+ int r = 0;
+ if(comm.Get_rank() == 0) {
+ r = multi_sim_comm.Get_rank();
+ }
+ comm.Bcast(r, 0);
+ if(r>0) {
+ gridfname_tmp = "/dev/null";
+ }
+ ofile->enforceSuffix("");
+ }
+ ofile->open(gridfname_tmp);
+ ofile->setHeavyFlush();
+ gridfiles_.push_back(ofile);
+ }
}
- ofile->setHeavyFlush();
- // output periodicities of variables
- ofile->setupPrintValue( getPntrToArgument(i) );
- // push back
- hillsOfiles_.push_back(ofile);
- }
log<<" Bibliography "<<plumed.cite("Pfaendtner and Bonomi. J. Chem. Theory Comput. 11, 5062 (2015)");
if(doInt_) log<<plumed.cite(
"Baftizadeh, Cossio, Pietrucci, and Laio, Curr. Phys. Chem. 2, 79 (2012)");
if(multiple_w) log<<plumed.cite(
"Raiteri, Laio, Gervasio, Micheletti, and Parrinello, J. Phys. Chem. B 110, 3533 (2006)");
-
log<<"\n";
-
- turnOnDerivatives();
}
void PBMetaD::readGaussians(int iarg, IFile *ifile){
@@ -619,6 +698,7 @@ double PBMetaD::getBiasAndDerivatives(int iarg, const vector<double>& cv, double
bias = BiasGrids_[iarg]->getValue(cv);
}
}
+
return bias;
}
@@ -672,7 +752,6 @@ void PBMetaD::calculate()
for(unsigned i=0; i<getNumberOfArguments(); ++i){
const double f = - exp(-bias[i]/kbt_) / (ene) * deriv[i];
setOutputForce(i, f);
- getPntrToComponent("bias")->addDerivative(i,-f);
}
delete [] der;
@@ -688,11 +767,9 @@ void PBMetaD::update()
if(getStep()%stride_==0 && !isFirstStep ){nowAddAHill=true;}else{nowAddAHill=false;isFirstStep=false;}
if(nowAddAHill){
-
// get all CVs value
vector<double> cv(getNumberOfArguments());
for(unsigned i=0; i<getNumberOfArguments(); ++i) cv[i] = getArgument(i);
-
// get all biases and heights
vector<double> bias(getNumberOfArguments());
vector<double> height(getNumberOfArguments());
@@ -711,7 +788,7 @@ void PBMetaD::update()
height[i] *= height0_ / norm;
if(welltemp_) height[i] *= exp(-bias[i]/(kbt_*(biasf_-1.0)));
}
-
+
// Multiple walkers: share hills and add them all
if(multiple_w){
int nw = 0;
@@ -760,7 +837,17 @@ void PBMetaD::update()
writeGaussian(i, newhill, hillsOfiles_[i]);
}
}
- }
+
+ // write grid files
+ if(wgridstride_>0 && getStep()%wgridstride_==0) {
+ for(unsigned i=0; i<gridfiles_.size(); ++i) {
+ gridfiles_[i]->rewind();
+ BiasGrids_[i]->writeToFile(*gridfiles_[i]);
+ gridfiles_[i]->flush();
+ }
+ }
+
+ }
}
/// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height
diff --git a/vim/vimsyntax.sh b/vim/vimsyntax.sh
index 4b6d0b4922f56775701cf8d324e3027e8bc5703b..b0d04cf26bf63d2c4b68a742bcaa30bfb8d0501b 100755
--- a/vim/vimsyntax.sh
+++ b/vim/vimsyntax.sh
@@ -26,8 +26,8 @@ endif
let b:current_syntax="plumedf"
if exists("g:plumed_shortcuts")
- map <F3> :PMinus<CR>
- map <F4> :PPlus<CR>
+ noremap <buffer> <F3> :PMinus<CR>
+ noremap <buffer> <F4> :PPlus<CR>
endif
command! -nargs=0 -count=1 PPlus call PlumedMove(<count>)
@@ -145,15 +145,15 @@ endif
let b:current_syntax="plumed"
if exists("g:plumed_shortcuts")
- noremap <F2> <Esc>:PHelp<CR>
- inoremap <F2> <C-o>:PHelp<CR>
+ noremap <buffer> <F2> :PHelp<CR>
+ inoremap <buffer> <F2> <Esc>:PHelp<CR>
endif
" path for plumed plugin
let s:path=expand('<sfile>:p:h:h')
" All except space and hash are in word
-set iskeyword=33,34,36-126
+setlocal iskeyword=33,34,36-126
" Matching dots, possibly followed by a comment
" Only dots are part of the match
@@ -329,6 +329,10 @@ fun! PlumedContextManual()
execute 'rightbelow split | view ' name
endif
let b:plumed_helpfile=1
+" this is to allow closing the window with F2
+ if exists("g:plumed_shortcuts")
+ noremap <buffer> <F2> :PHelp<CR>
+ endif
endif
endfunction