From a8cfda16466b0ab97ab73aa5cc0fc2586e3dd2c9 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sun, 1 May 2016 11:53:17 +0100
Subject: [PATCH] Started adding all the stuff for reweighting in an averaging
base class
---
.../rt-histo/analysis.0.fesA1.reference | 2 +-
regtest/analysis/rt-histo/fesA1.reference | 2 +-
regtest/analysis/rt-histo/fesA1u.reference | 2 +-
regtest/analysis/rt-histo/fesA2.reference | 2 +-
regtest/analysis/rt-histo/fesA2u.reference | 2 +-
regtest/analysis/rt-histo/fesA3.reference | 2 +-
regtest/analysis/rt0/histoA.reference | 20 +--
regtest/analysis/rt0/histoB.reference | 28 ++--
regtest/analysis/rt0/histoD.reference | 4 +-
regtest/analysis/rt0/histoD1.reference | 4 +-
regtest/analysis/rt0/histoE.reference | 4 +-
regtest/analysis/rt0/plumed.dat | 34 +++--
.../analysis.0.contour2.dat.reference | 2 +-
.../rt-wcsurface2/contour2.dat.reference | 2 +-
.../analysis.0.dens1b.reference | 2 +-
.../rt-dens-average/dens1.reference | 2 +-
.../rt-dens-average/dens1b.reference | 2 +-
.../rt-dens-average/dens1ba.reference | 104 +++++++--------
.../rt-dens-average/dens2.reference | 2 +-
.../rt-dens-average/dens5.reference | 2 +-
.../rt-dens-average/dens5a.reference | 104 +++++++--------
.../rt-dens-average/interpol.reference | 2 +-
.../multicolvar/rt-dens-average/plumed.dat | 20 ++-
.../multicolvar/rt-dens/interpol.reference | 2 +-
regtest/multicolvar/rt-dens/plumed.dat | 4 +-
src/analysis/Histogram.cpp | 22 +--
src/bias/ReweightBase.cpp | 73 ++++++++++
src/bias/ReweightBase.h | 51 +++++++
src/bias/ReweightBias.cpp | 67 ++++++++++
src/bias/ReweightMetad.cpp | 67 ++++++++++
src/bias/ReweightTemperature.cpp | 78 +++++++++++
src/core/ActionWithArguments.h | 2 +-
src/gridtools/ActionWithGrid.cpp | 88 ++----------
src/gridtools/ActionWithGrid.h | 33 +----
src/gridtools/ActionWithInputGrid.cpp | 17 ++-
src/gridtools/ActionWithInputGrid.h | 6 +-
src/gridtools/AverageOnGrid.cpp | 17 +--
src/gridtools/AverageOnGrid.h | 1 -
src/gridtools/ContourFindingBase.cpp | 2 +-
src/gridtools/ContourFindingBase.h | 2 +-
src/gridtools/ConvertToFES.cpp | 5 +-
src/gridtools/DumpCube.cpp | 42 +++---
src/gridtools/DumpGrid.cpp | 53 ++------
src/gridtools/FindContour.cpp | 12 +-
src/gridtools/FindContourSurface.cpp | 22 +--
src/gridtools/FourierTransform.cpp | 6 +-
src/gridtools/GridPrintingBase.cpp | 76 +++++++++++
src/gridtools/GridPrintingBase.h | 48 +++++++
src/gridtools/GridVessel.cpp | 46 ++-----
src/gridtools/GridVessel.h | 45 +------
src/gridtools/HistogramOnGrid.cpp | 13 +-
src/gridtools/HistogramOnGrid.h | 2 +-
src/gridtools/InterpolateGrid.cpp | 17 +--
src/multicolvar/MultiColvarDensity.cpp | 16 +--
src/vesselbase/ActionWithAveraging.cpp | 126 ++++++++++++++++++
src/vesselbase/ActionWithAveraging.h | 94 +++++++++++++
src/vesselbase/ActionWithVessel.cpp | 2 -
src/vesselbase/AveragingVessel.cpp | 62 +++++++++
src/vesselbase/AveragingVessel.h | 101 ++++++++++++++
59 files changed, 1157 insertions(+), 513 deletions(-)
create mode 100644 src/bias/ReweightBase.cpp
create mode 100644 src/bias/ReweightBase.h
create mode 100644 src/bias/ReweightBias.cpp
create mode 100644 src/bias/ReweightMetad.cpp
create mode 100644 src/bias/ReweightTemperature.cpp
create mode 100644 src/gridtools/GridPrintingBase.cpp
create mode 100644 src/gridtools/GridPrintingBase.h
create mode 100644 src/vesselbase/ActionWithAveraging.cpp
create mode 100644 src/vesselbase/ActionWithAveraging.h
create mode 100644 src/vesselbase/AveragingVessel.cpp
create mode 100644 src/vesselbase/AveragingVessel.h
diff --git a/regtest/analysis/rt-histo/analysis.0.fesA1.reference b/regtest/analysis/rt-histo/analysis.0.fesA1.reference
index 238a9b202..ba54f4733 100644
--- a/regtest/analysis/rt-histo/analysis.0.fesA1.reference
+++ b/regtest/analysis/rt-histo/analysis.0.fesA1.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1 dfA1_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt-histo/fesA1.reference b/regtest/analysis/rt-histo/fesA1.reference
index 6edf125bb..5b9a304f3 100644
--- a/regtest/analysis/rt-histo/fesA1.reference
+++ b/regtest/analysis/rt-histo/fesA1.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1 dfA1_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt-histo/fesA1u.reference b/regtest/analysis/rt-histo/fesA1u.reference
index c15e81f17..674191daf 100644
--- a/regtest/analysis/rt-histo/fesA1u.reference
+++ b/regtest/analysis/rt-histo/fesA1u.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1u dfA1u_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt-histo/fesA2.reference b/regtest/analysis/rt-histo/fesA2.reference
index 6edf125bb..5b9a304f3 100644
--- a/regtest/analysis/rt-histo/fesA2.reference
+++ b/regtest/analysis/rt-histo/fesA2.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1 dfA1_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt-histo/fesA2u.reference b/regtest/analysis/rt-histo/fesA2u.reference
index c15e81f17..674191daf 100644
--- a/regtest/analysis/rt-histo/fesA2u.reference
+++ b/regtest/analysis/rt-histo/fesA2u.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1u dfA1u_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt-histo/fesA3.reference b/regtest/analysis/rt-histo/fesA3.reference
index 6edf125bb..5b9a304f3 100644
--- a/regtest/analysis/rt-histo/fesA3.reference
+++ b/regtest/analysis/rt-histo/fesA3.reference
@@ -1,5 +1,5 @@
#! FIELDS x fA1 dfA1_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
diff --git a/regtest/analysis/rt0/histoA.reference b/regtest/analysis/rt0/histoA.reference
index a9c66815b..da854d9bc 100644
--- a/regtest/analysis/rt0/histoA.reference
+++ b/regtest/analysis/rt0/histoA.reference
@@ -44,28 +44,28 @@
1.1100 0.0000 0.0000
1.1400 0.0031 0.1101
1.1700 0.0086 0.2842
- 1.2000 0.0221 0.6648
+ 1.2000 0.0222 0.6648
1.2300 0.0521 1.4068
- 1.2600 0.1119 2.6873
+ 1.2600 0.1120 2.6873
1.2900 0.2199 4.6183
- 1.3200 0.3947 7.1055
+ 1.3200 0.3948 7.1055
1.3500 0.6476 9.7138
- 1.3800 0.9709 11.6511
- 1.4100 1.3304 11.9739
- 1.4400 1.6661 9.9968
+ 1.3800 0.9709 11.6512
+ 1.4100 1.3304 11.9738
+ 1.4400 1.6661 9.9967
1.4700 1.9069 5.7208
1.5000 1.9947 0.0000
1.5300 1.9069 -5.7208
- 1.5600 1.6661 -9.9968
- 1.5900 1.3304 -11.9739
+ 1.5600 1.6661 -9.9967
+ 1.5900 1.3304 -11.9738
1.6200 0.9724 -11.5957
1.6500 0.6520 -9.5611
- 1.6800 0.4067 -6.7240
+ 1.6800 0.4067 -6.7241
1.7100 0.2497 -3.7552
1.7400 0.1799 -0.9216
1.7700 0.1937 1.8508
1.8000 0.2921 4.7343
- 1.8300 0.4788 7.7100
+ 1.8300 0.4789 7.7099
1.8600 0.7517 10.3708
1.8900 1.0893 11.9819
1.9200 1.4485 11.5877
diff --git a/regtest/analysis/rt0/histoB.reference b/regtest/analysis/rt0/histoB.reference
index 5760efcba..7b569d51f 100644
--- a/regtest/analysis/rt0/histoB.reference
+++ b/regtest/analysis/rt0/histoB.reference
@@ -1,5 +1,5 @@
#! FIELDS x hB dhB_x
-#! SET normalisation 2.2220
+#! SET normalisation 4.0542
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
@@ -42,37 +42,37 @@
1.0500 0.0000 0.0000
1.0800 0.0000 0.0000
1.1100 0.0000 0.0000
- 1.1400 0.0027 0.0991
- 1.1700 0.0077 0.2558
+ 1.1400 0.0028 0.0991
+ 1.1700 0.0078 0.2558
1.2000 0.0199 0.5984
1.2300 0.0469 1.2663
1.2600 0.1008 2.4189
- 1.2900 0.1979 4.1570
+ 1.2900 0.1979 4.1569
1.3200 0.3553 6.3957
1.3500 0.5829 8.7435
- 1.3800 0.8740 10.4873
+ 1.3800 0.8739 10.4873
1.4100 1.1975 10.7777
1.4400 1.4997 8.9981
- 1.4700 1.7165 5.1493
- 1.5000 1.7954 0.0000
- 1.5300 1.7165 -5.1493
+ 1.4700 1.7164 5.1493
+ 1.5000 1.7955 0.0000
+ 1.5300 1.7164 -5.1493
1.5600 1.4997 -8.9981
1.5900 1.1975 -10.7777
- 1.6200 0.8756 -10.4263
+ 1.6200 0.8755 -10.4263
1.6500 0.5877 -8.5755
1.6800 0.3684 -5.9762
1.7100 0.2307 -3.2076
1.7400 0.1755 -0.4767
1.7700 0.2027 2.3167
- 1.8000 0.3169 5.3400
+ 1.8000 0.3169 5.3401
1.8300 0.5250 8.5369
- 1.8600 0.8262 11.4287
+ 1.8600 0.8262 11.4288
1.8900 1.1981 13.1788
1.9200 1.5931 12.7452
- 1.9500 1.9362 9.6808
- 1.9800 2.1505 4.3010
+ 1.9500 1.9362 9.6809
+ 1.9800 2.1505 4.3011
2.0100 2.1830 -2.1830
- 2.0400 2.0253 -8.1011
+ 2.0400 2.0253 -8.1012
2.0700 1.7172 -12.0206
2.1000 1.3307 -13.3071
2.1300 0.9424 -12.2517
diff --git a/regtest/analysis/rt0/histoD.reference b/regtest/analysis/rt0/histoD.reference
index c2827416c..a05fbfe36 100644
--- a/regtest/analysis/rt0/histoD.reference
+++ b/regtest/analysis/rt0/histoD.reference
@@ -1,5 +1,5 @@
#! FIELDS x hD
-#! SET normalisation 5.5474
+#! SET normalisation 4.0542
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
@@ -37,7 +37,7 @@
0.9000 0.0000
0.9300 0.0000
0.9600 0.0000
- 0.9900 1.4932
+ 0.9900 0.0000
1.0200 0.0000
1.0500 0.0000
1.0800 0.0000
diff --git a/regtest/analysis/rt0/histoD1.reference b/regtest/analysis/rt0/histoD1.reference
index 7a6242f4e..10fd0663e 100644
--- a/regtest/analysis/rt0/histoD1.reference
+++ b/regtest/analysis/rt0/histoD1.reference
@@ -1,5 +1,5 @@
#! FIELDS x hD1
-#! SET normalisation 1.4932
+#! SET normalisation 0.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
@@ -37,7 +37,7 @@
0.9000 0.0000
0.9300 0.0000
0.9600 0.0000
- 0.9900 1.4932
+ 0.9900 0.0000
1.0200 0.0000
1.0500 0.0000
1.0800 0.0000
diff --git a/regtest/analysis/rt0/histoE.reference b/regtest/analysis/rt0/histoE.reference
index b36b2558f..1f18687de 100644
--- a/regtest/analysis/rt0/histoE.reference
+++ b/regtest/analysis/rt0/histoE.reference
@@ -1,5 +1,5 @@
#! FIELDS x fE
-#! SET normalisation 3.0547
+#! SET normalisation 1.0000
#! SET min_x 0.0
#! SET max_x 3.0
#! SET nbins_x 100
@@ -37,7 +37,7 @@
0.9000 inf
0.9300 inf
0.9600 inf
- 0.9900 -1.0000
+ 0.9900 inf
1.0200 inf
1.0500 inf
1.0800 inf
diff --git a/regtest/analysis/rt0/plumed.dat b/regtest/analysis/rt0/plumed.dat
index b5ff3b6cd..db0eeb44f 100755
--- a/regtest/analysis/rt0/plumed.dat
+++ b/regtest/analysis/rt0/plumed.dat
@@ -4,7 +4,6 @@ RESTRAINT ARG=x SLOPE=1.0 AT=0.0
HISTOGRAM ...
ARG=x
- TEMP=300
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
@@ -12,20 +11,21 @@ HISTOGRAM ...
LABEL=hA
... HISTOGRAM
-PRINT_GRID GRID=hA FILE=histoA STRIDE=1 FMT=%8.4f
+DUMPGRID GRID=hA FILE=histoA STRIDE=1 FMT=%8.4f
+
+bias: REWEIGHT_BIAS TEMP=300
HISTOGRAM ...
ARG=x
- TEMP=300
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
BANDWIDTH=0.1
- REWEIGHT_BIAS
+ LOGWEIGHTS=bias
LABEL=hB
... HISTOGRAM
-PRINT_GRID GRID=hB FILE=histoB STRIDE=1 FMT=%8.4f
+DUMPGRID GRID=hB FILE=histoB STRIDE=1 FMT=%8.4f
HISTOGRAM ...
ARG=x
@@ -36,68 +36,66 @@ HISTOGRAM ...
LABEL=hC
... HISTOGRAM
-PRINT_GRID GRID=hC FILE=histoC USE_ALL_DATA
+DUMPGRID GRID=hC FILE=histoC
HISTOGRAM ...
ARG=x
- TEMP=300
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
KERNEL=DISCRETE
- REWEIGHT_BIAS
+ LOGWEIGHTS=bias
UNORMALIZED
LABEL=hD
... HISTOGRAM
-PRINT_GRID GRID=hD FILE=histoD FMT=%8.4f USE_ALL_DATA
+DUMPGRID GRID=hD FILE=histoD FMT=%8.4f
+
+bias-ht: REWEIGHT_BIAS TEMP=10000
HISTOGRAM ...
ARG=x
- TEMP=10000
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
KERNEL=DISCRETE
- REWEIGHT_BIAS
+ LOGWEIGHTS=bias-ht
UNORMALIZED
LABEL=hE
... HISTOGRAM
fE: CONVERT_TO_FES GRID=hE TEMP=10000
-PRINT_GRID GRID=fE FMT=%8.4f USE_ALL_DATA FILE=histoE
+DUMPGRID GRID=fE FMT=%8.4f FILE=histoE
HISTOGRAM ...
ARG=x
- TEMP=300
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
KERNEL=DISCRETE
- REWEIGHT_BIAS
+ LOGWEIGHTS=bias
UPDATE_FROM=0
UPDATE_UNTIL=1
UNORMALIZED
LABEL=hD1
... HISTOGRAM
-PRINT_GRID GRID=hD1 FILE=histoD1 FMT=%8.4f USE_ALL_DATA
+DUMPGRID GRID=hD1 FILE=histoD1 FMT=%8.4f
HISTOGRAM ...
ARG=x
- TEMP=300
GRID_MIN=0.0
GRID_MAX=3.0
GRID_BIN=100
KERNEL=DISCRETE
- REWEIGHT_BIAS
+ LOGWEIGHTS=bias
UPDATE_FROM=1
UPDATE_UNTIL=2
UNORMALIZED
LABEL=hD2
... HISTOGRAM
-PRINT_GRID GRID=hD2 FILE=histoD2 FMT=%8.4f USE_ALL_DATA
+DUMPGRID GRID=hD2 FILE=histoD2 FMT=%8.4f
diff --git a/regtest/crystallization/rt-wcsurface2/analysis.0.contour2.dat.reference b/regtest/crystallization/rt-wcsurface2/analysis.0.contour2.dat.reference
index bdc34b0f0..624036870 100644
--- a/regtest/crystallization/rt-wcsurface2/analysis.0.contour2.dat.reference
+++ b/regtest/crystallization/rt-wcsurface2/analysis.0.contour2.dat.reference
@@ -1,5 +1,5 @@
#! FIELDS x y ss2
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -7.27971
#! SET max_x 7.27971
#! SET nbins_x 14
diff --git a/regtest/crystallization/rt-wcsurface2/contour2.dat.reference b/regtest/crystallization/rt-wcsurface2/contour2.dat.reference
index bdc34b0f0..624036870 100644
--- a/regtest/crystallization/rt-wcsurface2/contour2.dat.reference
+++ b/regtest/crystallization/rt-wcsurface2/contour2.dat.reference
@@ -1,5 +1,5 @@
#! FIELDS x y ss2
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -7.27971
#! SET max_x 7.27971
#! SET nbins_x 14
diff --git a/regtest/multicolvar/rt-dens-average/analysis.0.dens1b.reference b/regtest/multicolvar/rt-dens-average/analysis.0.dens1b.reference
index 7a16efe95..f1a1ba1c4 100644
--- a/regtest/multicolvar/rt-dens-average/analysis.0.dens1b.reference
+++ b/regtest/multicolvar/rt-dens-average/analysis.0.dens1b.reference
@@ -1,5 +1,5 @@
#! FIELDS x dens.dens ddens.dens_x density ddensity_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
diff --git a/regtest/multicolvar/rt-dens-average/dens1.reference b/regtest/multicolvar/rt-dens-average/dens1.reference
index 147b4e1b3..0f0bb8187 100644
--- a/regtest/multicolvar/rt-dens-average/dens1.reference
+++ b/regtest/multicolvar/rt-dens-average/dens1.reference
@@ -1,5 +1,5 @@
#! FIELDS x dens.dens ddens.dens_x density ddensity_x
-#! SET normalisation 0.0000
+#! SET normalisation 2.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
diff --git a/regtest/multicolvar/rt-dens-average/dens1b.reference b/regtest/multicolvar/rt-dens-average/dens1b.reference
index 9fb24891b..074db3f29 100644
--- a/regtest/multicolvar/rt-dens-average/dens1b.reference
+++ b/regtest/multicolvar/rt-dens-average/dens1b.reference
@@ -1,5 +1,5 @@
#! FIELDS x dens.dens ddens.dens_x density ddensity_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
diff --git a/regtest/multicolvar/rt-dens-average/dens1ba.reference b/regtest/multicolvar/rt-dens-average/dens1ba.reference
index a92621e9d..d60f11274 100644
--- a/regtest/multicolvar/rt-dens-average/dens1ba.reference
+++ b/regtest/multicolvar/rt-dens-average/dens1ba.reference
@@ -1,56 +1,56 @@
-#! FIELDS x dens.dens ddens.dens_x
-#! SET normalisation 0.0000
+#! FIELDS x dens.dens ddens.dens_x density ddensity_x
+#! SET normalisation 1.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
#! SET periodic_x true
- -2.5000 1.0000 0.0000
- -2.4000 1.0000 0.0000
- -2.3000 1.0000 0.0000
- -2.2000 1.0000 0.0000
- -2.1000 0.0000 0.0000
- -2.0000 0.0000 0.0000
- -1.9000 0.0000 0.0000
- -1.8000 0.0000 0.0000
- -1.7000 0.0000 0.0000
- -1.6000 0.0000 0.0000
- -1.5000 0.0000 0.0000
- -1.4000 0.0000 0.0000
- -1.3000 0.0000 0.0000
- -1.2000 0.0000 0.0000
- -1.1000 0.0000 0.0000
- -1.0000 0.0000 0.0000
- -0.9000 0.0000 0.0000
- -0.8000 0.0000 0.0000
- -0.7000 0.0000 0.0000
- -0.6000 0.0000 0.0000
- -0.5000 0.0000 0.0000
- -0.4000 0.0000 0.0000
- -0.3000 0.0000 0.0000
- -0.2000 0.0000 0.0000
- -0.1000 0.0000 0.0000
- 0.0000 0.0000 0.0000
- 0.1000 0.0000 0.0000
- 0.2000 0.0000 0.0000
- 0.3000 0.0000 0.0000
- 0.4000 0.0000 0.0000
- 0.5000 0.0000 0.0000
- 0.6000 0.0000 0.0000
- 0.7000 0.0000 0.0000
- 0.8000 0.0000 0.0000
- 0.9000 0.0000 0.0000
- 1.0000 0.0000 0.0000
- 1.1000 0.0000 0.0000
- 1.2000 1.0000 0.0000
- 1.3000 1.0000 0.0000
- 1.4000 1.0000 0.0000
- 1.5000 1.0000 -0.0000
- 1.6000 1.0000 0.0000
- 1.7000 1.0000 0.0000
- 1.8000 1.0000 0.0000
- 1.9000 1.0000 0.0000
- 2.0000 1.0000 0.0000
- 2.1000 1.0000 0.0000
- 2.2000 1.0000 -0.0000
- 2.3000 1.0000 0.0000
- 2.4000 1.0000 -0.0000
+ -2.5000 1.0000 0.0000 0.0876 -1.0955
+ -2.4000 1.0000 0.0000 0.0222 -0.3324
+ -2.3000 1.0000 0.0000 0.0044 -0.0764
+ -2.2000 1.0000 0.0000 0.0007 -0.0134
+ -2.1000 0.0000 0.0000 0.0000 0.0000
+ -2.0000 0.0000 0.0000 0.0000 0.0000
+ -1.9000 0.0000 0.0000 0.0000 0.0000
+ -1.8000 0.0000 0.0000 0.0000 0.0000
+ -1.7000 0.0000 0.0000 0.0000 0.0000
+ -1.6000 0.0000 0.0000 0.0000 0.0000
+ -1.5000 0.0000 0.0000 0.0000 0.0000
+ -1.4000 0.0000 0.0000 0.0000 0.0000
+ -1.3000 0.0000 0.0000 0.0000 0.0000
+ -1.2000 0.0000 0.0000 0.0000 0.0000
+ -1.1000 0.0000 0.0000 0.0000 0.0000
+ -1.0000 0.0000 0.0000 0.0000 0.0000
+ -0.9000 0.0000 0.0000 0.0000 0.0000
+ -0.8000 0.0000 0.0000 0.0000 0.0000
+ -0.7000 0.0000 0.0000 0.0000 0.0000
+ -0.6000 0.0000 0.0000 0.0000 0.0000
+ -0.5000 0.0000 0.0000 0.0000 0.0000
+ -0.4000 0.0000 0.0000 0.0000 0.0000
+ -0.3000 0.0000 0.0000 0.0000 0.0000
+ -0.2000 0.0000 0.0000 0.0000 0.0000
+ -0.1000 0.0000 0.0000 0.0000 0.0000
+ 0.0000 0.0000 0.0000 0.0000 0.0000
+ 0.1000 0.0000 0.0000 0.0000 0.0000
+ 0.2000 0.0000 0.0000 0.0000 0.0000
+ 0.3000 0.0000 0.0000 0.0000 0.0000
+ 0.4000 0.0000 0.0000 0.0000 0.0000
+ 0.5000 0.0000 0.0000 0.0000 0.0000
+ 0.6000 0.0000 0.0000 0.0000 0.0000
+ 0.7000 0.0000 0.0000 0.0000 0.0000
+ 0.8000 0.0000 0.0000 0.0000 0.0000
+ 0.9000 0.0000 0.0000 0.0000 0.0000
+ 1.0000 0.0000 0.0000 0.0000 0.0000
+ 1.1000 0.0000 0.0000 0.0000 0.0000
+ 1.2000 1.0000 0.0000 0.0007 0.0134
+ 1.3000 1.0000 0.0000 0.0044 0.0764
+ 1.4000 1.0000 0.0000 0.0222 0.3324
+ 1.5000 1.0000 -0.0000 0.0876 1.0955
+ 1.6000 1.0000 0.0000 0.2700 2.6995
+ 1.7000 1.0000 0.0000 0.6476 4.8569
+ 1.8000 1.0000 0.0000 1.2099 6.0493
+ 1.9000 1.0000 0.0000 1.7603 4.4008
+ 2.0000 1.0000 0.0000 1.9947 0.0000
+ 2.1000 1.0000 0.0000 1.7603 -4.4008
+ 2.2000 1.0000 -0.0000 1.2099 -6.0493
+ 2.3000 1.0000 0.0000 0.6476 -4.8569
+ 2.4000 1.0000 -0.0000 0.2700 -2.6995
diff --git a/regtest/multicolvar/rt-dens-average/dens2.reference b/regtest/multicolvar/rt-dens-average/dens2.reference
index 147b4e1b3..0f0bb8187 100644
--- a/regtest/multicolvar/rt-dens-average/dens2.reference
+++ b/regtest/multicolvar/rt-dens-average/dens2.reference
@@ -1,5 +1,5 @@
#! FIELDS x dens.dens ddens.dens_x density ddensity_x
-#! SET normalisation 0.0000
+#! SET normalisation 2.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
diff --git a/regtest/multicolvar/rt-dens-average/dens5.reference b/regtest/multicolvar/rt-dens-average/dens5.reference
index 147b4e1b3..0f0bb8187 100644
--- a/regtest/multicolvar/rt-dens-average/dens5.reference
+++ b/regtest/multicolvar/rt-dens-average/dens5.reference
@@ -1,5 +1,5 @@
#! FIELDS x dens.dens ddens.dens_x density ddensity_x
-#! SET normalisation 0.0000
+#! SET normalisation 2.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
diff --git a/regtest/multicolvar/rt-dens-average/dens5a.reference b/regtest/multicolvar/rt-dens-average/dens5a.reference
index 5cd1d702e..a34c68c87 100644
--- a/regtest/multicolvar/rt-dens-average/dens5a.reference
+++ b/regtest/multicolvar/rt-dens-average/dens5a.reference
@@ -1,56 +1,56 @@
-#! FIELDS x dens.dens ddens.dens_x
-#! SET normalisation 0.0000
+#! FIELDS x dens.dens ddens.dens_x density ddensity_x
+#! SET normalisation 2.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 50
#! SET periodic_x true
- -2.5000 1.0000 0.0000
- -2.4000 1.0000 0.0000
- -2.3000 1.0000 0.0000
- -2.2000 1.0000 0.0000
- -2.1000 0.0000 0.0000
- -2.0000 0.0000 0.0000
- -1.9000 0.0000 0.0000
- -1.8000 0.0000 0.0000
- -1.7000 0.0000 0.0000
- -1.6000 0.0000 0.0000
- -1.5000 0.0000 0.0000
- -1.4000 0.0000 0.0000
- -1.3000 0.0000 0.0000
- -1.2000 0.0000 0.0000
- -1.1000 0.0000 0.0000
- -1.0000 0.0000 0.0000
- -0.9000 0.0000 0.0000
- -0.8000 0.0000 0.0000
- -0.7000 0.0000 0.0000
- -0.6000 0.0000 0.0000
- -0.5000 0.0000 0.0000
- -0.4000 0.0000 0.0000
- -0.3000 0.0000 0.0000
- -0.2000 0.0000 0.0000
- -0.1000 0.0000 0.0000
- 0.0000 0.0000 0.0000
- 0.1000 0.0000 0.0000
- 0.2000 2.0000 -0.0000
- 0.3000 2.0000 0.0000
- 0.4000 2.0000 0.0000
- 0.5000 2.0000 -0.0000
- 0.6000 2.0000 0.0000
- 0.7000 2.0000 0.0000
- 0.8000 2.0000 0.0000
- 0.9000 2.0000 -0.0000
- 1.0000 2.0000 0.0000
- 1.1000 2.0000 0.0000
- 1.2000 1.9994 -0.0138
- 1.3000 1.9933 -0.1662
- 1.4000 1.9241 -1.7526
- 1.5000 1.5000 -6.2500
- 1.6000 1.0759 -1.7526
- 1.7000 1.0067 -0.1662
- 1.8000 1.0006 -0.0138
- 1.9000 1.0000 0.0000
- 2.0000 1.0000 0.0000
- 2.1000 1.0000 0.0000
- 2.2000 1.0000 -0.0000
- 2.3000 1.0000 0.0000
- 2.4000 1.0000 -0.0000
+ -2.5000 1.0000 0.0000 0.0438 -0.5478
+ -2.4000 1.0000 0.0000 0.0111 -0.1662
+ -2.3000 1.0000 0.0000 0.0022 -0.0382
+ -2.2000 1.0000 0.0000 0.0003 -0.0067
+ -2.1000 0.0000 0.0000 0.0000 0.0000
+ -2.0000 0.0000 0.0000 0.0000 0.0000
+ -1.9000 0.0000 0.0000 0.0000 0.0000
+ -1.8000 0.0000 0.0000 0.0000 0.0000
+ -1.7000 0.0000 0.0000 0.0000 0.0000
+ -1.6000 0.0000 0.0000 0.0000 0.0000
+ -1.5000 0.0000 0.0000 0.0000 0.0000
+ -1.4000 0.0000 0.0000 0.0000 0.0000
+ -1.3000 0.0000 0.0000 0.0000 0.0000
+ -1.2000 0.0000 0.0000 0.0000 0.0000
+ -1.1000 0.0000 0.0000 0.0000 0.0000
+ -1.0000 0.0000 0.0000 0.0000 0.0000
+ -0.9000 0.0000 0.0000 0.0000 0.0000
+ -0.8000 0.0000 0.0000 0.0000 0.0000
+ -0.7000 0.0000 0.0000 0.0000 0.0000
+ -0.6000 0.0000 0.0000 0.0000 0.0000
+ -0.5000 0.0000 0.0000 0.0000 0.0000
+ -0.4000 0.0000 0.0000 0.0000 0.0000
+ -0.3000 0.0000 0.0000 0.0000 0.0000
+ -0.2000 0.0000 0.0000 0.0000 0.0000
+ -0.1000 0.0000 0.0000 0.0000 0.0000
+ 0.0000 0.0000 0.0000 0.0000 0.0000
+ 0.1000 0.0000 0.0000 0.0000 0.0000
+ 0.2000 2.0000 -0.0000 0.0003 0.0067
+ 0.3000 2.0000 0.0000 0.0022 0.0382
+ 0.4000 2.0000 0.0000 0.0111 0.1662
+ 0.5000 2.0000 -0.0000 0.0438 0.5478
+ 0.6000 2.0000 0.0000 0.1350 1.3498
+ 0.7000 2.0000 0.0000 0.3238 2.4285
+ 0.8000 2.0000 0.0000 0.6049 3.0246
+ 0.9000 2.0000 -0.0000 0.8802 2.2004
+ 1.0000 2.0000 0.0000 0.9974 0.0000
+ 1.1000 2.0000 0.0000 0.8802 -2.2004
+ 1.2000 1.9994 -0.0138 0.6053 -3.0179
+ 1.3000 1.9933 -0.1662 0.3260 -2.3903
+ 1.4000 1.9241 -1.7526 0.1461 -1.1836
+ 1.5000 1.5000 -6.2500 0.0876 -0.0000
+ 1.6000 1.0759 -1.7526 0.1461 1.1836
+ 1.7000 1.0067 -0.1662 0.3260 2.3903
+ 1.8000 1.0006 -0.0138 0.6053 3.0179
+ 1.9000 1.0000 0.0000 0.8802 2.2004
+ 2.0000 1.0000 0.0000 0.9974 0.0000
+ 2.1000 1.0000 0.0000 0.8802 -2.2004
+ 2.2000 1.0000 -0.0000 0.6049 -3.0246
+ 2.3000 1.0000 0.0000 0.3238 -2.4285
+ 2.4000 1.0000 -0.0000 0.1350 -1.3498
diff --git a/regtest/multicolvar/rt-dens-average/interpol.reference b/regtest/multicolvar/rt-dens-average/interpol.reference
index 75d4ea449..3ff0707ad 100644
--- a/regtest/multicolvar/rt-dens-average/interpol.reference
+++ b/regtest/multicolvar/rt-dens-average/interpol.reference
@@ -1,5 +1,5 @@
#! FIELDS x interpol dinterpol_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 100
diff --git a/regtest/multicolvar/rt-dens-average/plumed.dat b/regtest/multicolvar/rt-dens-average/plumed.dat
index 371be5092..30bba5a3c 100644
--- a/regtest/multicolvar/rt-dens-average/plumed.dat
+++ b/regtest/multicolvar/rt-dens-average/plumed.dat
@@ -3,22 +3,28 @@ UNITS NATURAL
dens: COORDINATIONNUMBER SPECIESA=2 SPECIESB=3-5 SWITCH={RATIONAL D_0=1.1 R_0=0.001 D_MAX=2.0}
# Print the average (whole trajectory) density with a stride of one
-dens1: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 BANDWIDTH=0.2
+dens1: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 UNORMALIZED BANDWIDTH=0.2
DUMPGRID GRID=dens1 STRIDE=1 FILE=dens1 FMT=%8.4f
# Print the average density with a stride of two
DUMPGRID GRID=dens1 STRIDE=2 FILE=dens2 FMT=%8.4f
# Print the average density (whole trajectory)
DUMPGRID GRID=dens1 FILE=dens5 FMT=%8.4f
-DUMPGRID GRID=dens1 FILE=dens5a PRINT_AVERAGE FMT=%8.4f
-# Print block averages (over two frames) of the density
-dens1b: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 CLEAR=1 BANDWIDTH=0.2
+# Print the average (whole trajectory) density with a stride of one
+dens1a: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 BANDWIDTH=0.2
+DUMPGRID GRID=dens1a FILE=dens5a FMT=%8.4f
+
+# Print unormalized averages (over two frames) of the density
+dens1b: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 CLEAR=1 UNORMALIZED BANDWIDTH=0.2
DUMPGRID GRID=dens1b STRIDE=1 FILE=dens1b FMT=%8.4f
-DUMPGRID GRID=dens1b STRIDE=1 FILE=dens1ba PRINT_AVERAGE FMT=%8.4f
+
+# Print block averages (over two frames) of the density
+dens1ba: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 CLEAR=1 BANDWIDTH=0.2
+DUMPGRID GRID=dens1ba STRIDE=1 FILE=dens1ba FMT=%8.4f
# Interpolate onto a finer grid and test
fine: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=100 BANDWIDTH=0.1
-interpol: INTERPOLATE_GRID GRID=dens1 GRID_BIN=100
+interpol: INTERPOLATE_GRID GRID=dens1a STRIDE=2 GRID_BIN=100
-DUMPGRID GRID=fine FILE=fine_grid PRINT_AVERAGE FMT=%8.4f
+DUMPGRID GRID=fine FILE=fine_grid FMT=%8.4f
DUMPGRID GRID=interpol FILE=interpol FMT=%8.4f
diff --git a/regtest/multicolvar/rt-dens/interpol.reference b/regtest/multicolvar/rt-dens/interpol.reference
index 851a78ece..70c0d6189 100644
--- a/regtest/multicolvar/rt-dens/interpol.reference
+++ b/regtest/multicolvar/rt-dens/interpol.reference
@@ -1,5 +1,5 @@
#! FIELDS x interpol dinterpol_x
-#! SET normalisation 0.0000
+#! SET normalisation 1.0000
#! SET min_x -2.5
#! SET max_x 2.5
#! SET nbins_x 100
diff --git a/regtest/multicolvar/rt-dens/plumed.dat b/regtest/multicolvar/rt-dens/plumed.dat
index b4094dd1e..828d9bb5d 100644
--- a/regtest/multicolvar/rt-dens/plumed.dat
+++ b/regtest/multicolvar/rt-dens/plumed.dat
@@ -13,8 +13,8 @@ dens2b: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=50 CLEAR=2 BANDWIDTH=0.1
DUMPGRID GRID=dens2b STRIDE=2 FILE=dens2b FMT=%8.4f
# Interpolate onto a finer grid and test
-fine: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=100 BANDWIDTH=0.1
-interpol: INTERPOLATE_GRID GRID=dens2 GRID_BIN=100
+fine: MULTICOLVARDENS DATA=dens ORIGIN=1 DIR=x NBINS=100 STRIDE=4 BANDWIDTH=0.1
+interpol: INTERPOLATE_GRID GRID=dens2 STRIDE=4 GRID_BIN=100
DUMPGRID GRID=fine FILE=fine_grid FMT=%8.4f
DUMPGRID GRID=interpol FILE=interpol FMT=%8.4f
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index f27c6535a..3a878a28b 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -117,8 +117,9 @@ private:
public:
static void registerKeywords( Keywords& keys );
explicit Histogram(const ActionOptions&ao);
- bool prepareForTasks();
- void finishTaskSet();
+ void prepareForAveraging();
+ void performOperations( const bool& from_update );
+ void finishAveraging();
bool isPeriodic(){ return false; }
unsigned getNumberOfDerivatives(){ return getNumberOfArguments(); }
void compute( const unsigned& , MultiValue& ) const ;
@@ -132,6 +133,7 @@ void Histogram::registerKeywords( Keywords& keys ){
keys.add("compulsory","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.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL");
}
Histogram::Histogram(const ActionOptions&ao):
@@ -166,31 +168,29 @@ kernel(NULL)
plumed_assert( myhist ); myhist->addOneKernelEachTimeOnly();
// Create a task list
for(unsigned i=0;i<mygrid->getNumberOfPoints();++i) addTaskToList(i);
- finishGridSetup(); resizeFunctions(); checkRead();
+ setAveragingAction( mygrid, myhist->noDiscreteKernels() ); checkRead();
}
-bool Histogram::prepareForTasks(){
+void Histogram::prepareForAveraging(){
// Now fetch the kernel and the active points
std::vector<double> point( getNumberOfArguments() );
for(unsigned i=0;i<point.size();++i) point[i]=getArgument(i);
- unsigned num_neigh; std::vector<unsigned> neighbors;
+ unsigned num_neigh; std::vector<unsigned> neighbors(1);
kernel = myhist->getKernelAndNeighbors( point, num_neigh, neighbors );
if( num_neigh>1 ){
// Activate relevant tasks
deactivateAllTasks();
for(unsigned i=0;i<num_neigh;++i) taskFlags[neighbors[i]]=1;
lockContributors();
-
- // Calculate the grid at all relevant points
- return true;
} else {
// This is used when we are not doing kernel density evaluation
- mygrid->setGridElement( neighbors[0], 0, mygrid->getGridElement( neighbors[0], 0 ) + 1.0 );
- return false;
+ mygrid->setGridElement( neighbors[0], 0, mygrid->getGridElement( neighbors[0], 0 ) + cweight );
}
}
-void Histogram::finishTaskSet(){
+void Histogram::performOperations( const bool& from_update ){ plumed_dbg_assert( !myhist->noDiscreteKernels() ); }
+
+void Histogram::finishAveraging(){
delete kernel;
}
diff --git a/src/bias/ReweightBase.cpp b/src/bias/ReweightBase.cpp
new file mode 100644
index 000000000..e0d13913e
--- /dev/null
+++ b/src/bias/ReweightBase.cpp
@@ -0,0 +1,73 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
+#include "core/Atoms.h"
+#include "ReweightBase.h"
+
+namespace PLMD {
+namespace bias {
+
+void ReweightBase::registerKeywords(Keywords& keys){
+ Action::registerKeywords( keys );
+ ActionWithValue::registerKeywords( keys );
+ ActionWithArguments::registerKeywords( keys );
+ keys.setComponentsIntroduction("This action calculates the logarithm of a weight for reweighting");
+ keys.add("optional","TEMP","the system temperature. This is not required if your MD code passes this quantity to PLUMED");
+ keys.remove("NUMERICAL_DERIVATIVES");
+}
+
+ReweightBase::ReweightBase(const ActionOptions&ao):
+Action(ao),
+ActionWithValue(ao),
+ActionWithArguments(ao)
+{
+ simtemp=0.; parse("TEMP",simtemp);
+ if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
+ else simtemp=plumed.getAtoms().getKbT();
+ if(simtemp==0) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP");
+ // Create something to hold the weight
+ addValue(); setNotPeriodic();
+}
+
+void ReweightBase::retrieveAllBiases( const std::string& lab, std::vector<Value*>& vals ){
+ std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
+ if( all.empty() ) error("your input file is not telling plumed to calculate anything");
+ log.printf(" using the following biases in reweighting ");
+ for(unsigned j=0;j<all.size();j++){
+ std::string flab; flab=all[j]->getLabel() + "." + lab;
+ if( all[j]->exists(flab) ){
+ vals.push_back( all[j]->copyOutput(flab) );
+ log.printf(" %s", flab.c_str());
+ }
+ }
+ log.printf("\n");
+ if( !vals.empty() ) requestArguments( vals );
+}
+
+void ReweightBase::calculate(){
+ double weight = getLogWeight();
+ setValue( weight );
+}
+
+}
+}
diff --git a/src/bias/ReweightBase.h b/src/bias/ReweightBase.h
new file mode 100644
index 000000000..ce982b37d
--- /dev/null
+++ b/src/bias/ReweightBase.h
@@ -0,0 +1,51 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_bias_ReweightBase_h
+#define __PLUMED_bias_ReweightBase_h
+
+#include "core/ActionWithValue.h"
+#include "core/ActionWithArguments.h"
+
+namespace PLMD {
+namespace bias {
+
+class ReweightBase :
+public ActionWithValue,
+public ActionWithArguments
+{
+protected:
+/// The temperature at which you are running the simulation
+ double simtemp;
+/// Retrieve all the values with a label ending in .lab
+ void retrieveAllBiases( const std::string& lab, std::vector<Value*>& vals );
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightBase(const ActionOptions&ao);
+ unsigned getNumberOfDerivatives(){ return 0; }
+ void calculate();
+ virtual double getLogWeight() const = 0;
+ void apply(){}
+};
+
+}
+}
+#endif
diff --git a/src/bias/ReweightBias.cpp b/src/bias/ReweightBias.cpp
new file mode 100644
index 000000000..db022301b
--- /dev/null
+++ b/src/bias/ReweightBias.cpp
@@ -0,0 +1,67 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "ReweightBase.h"
+
+//+PLUMEDOC REWEIGHTING REWEIGHT_BIAS
+/*
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace bias {
+
+class ReweightBias : public ReweightBase {
+private:
+/// The biases we are using in reweighting and the args we store them separately
+ std::vector<Value*> biases;
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightBias(const ActionOptions&ao);
+ double getLogWeight() const ;
+};
+
+PLUMED_REGISTER_ACTION(ReweightBias,"REWEIGHT_BIAS")
+
+void ReweightBias::registerKeywords(Keywords& keys ){
+ ReweightBase::registerKeywords( keys );
+}
+
+ReweightBias::ReweightBias(const ActionOptions&ao):
+Action(ao),
+ReweightBase(ao)
+{
+ retrieveAllBiases( "bias", biases );
+ if( biases.empty() ) error("there does not appear to be a bias acting on your system");
+}
+
+double ReweightBias::getLogWeight() const {
+ // Retrieve the bias
+ double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
+ return bias / simtemp;
+}
+
+}
+}
diff --git a/src/bias/ReweightMetad.cpp b/src/bias/ReweightMetad.cpp
new file mode 100644
index 000000000..78f66d52e
--- /dev/null
+++ b/src/bias/ReweightMetad.cpp
@@ -0,0 +1,67 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "ReweightBase.h"
+
+//+PLUMEDOC REWEIGHTING REWEIGHT_METAD
+/*
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace bias {
+
+class ReweightMetad : public ReweightBase {
+private:
+/// The biases we are using in reweighting and the args we store them separately
+ std::vector<Value*> biases;
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightMetad(const ActionOptions&ao);
+ double getLogWeight() const ;
+};
+
+PLUMED_REGISTER_ACTION(ReweightMetad,"REWEIGHT_METAD")
+
+void ReweightMetad::registerKeywords(Keywords& keys ){
+ ReweightBase::registerKeywords( keys );
+}
+
+ReweightMetad::ReweightMetad(const ActionOptions&ao):
+Action(ao),
+ReweightBase(ao)
+{
+ retrieveAllBiases( "rbias", biases );
+ if( biases.empty() ) error("there does not appear to be a metadynamics bias acting on your system");
+}
+
+double ReweightMetad::getLogWeight() const {
+ // Retrieve the bias
+ double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
+ return bias / simtemp;
+}
+
+}
+}
diff --git a/src/bias/ReweightTemperature.cpp b/src/bias/ReweightTemperature.cpp
new file mode 100644
index 000000000..e503ebfd7
--- /dev/null
+++ b/src/bias/ReweightTemperature.cpp
@@ -0,0 +1,78 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2011-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "core/ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/Atoms.h"
+#include "ReweightBase.h"
+
+//+PLUMEDOC REWEIGHTING REWEIGHT_TEMP
+/*
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
+namespace PLMD {
+namespace bias {
+
+class ReweightTemperature : public ReweightBase {
+private:
+///
+ double rtemp;
+/// The biases we are using in reweighting and the args we store them separately
+ std::vector<Value*> biases;
+public:
+ static void registerKeywords(Keywords&);
+ explicit ReweightTemperature(const ActionOptions&ao);
+ double getLogWeight() const ;
+};
+
+PLUMED_REGISTER_ACTION(ReweightTemperature,"REWEIGHT_TEMP")
+
+void ReweightTemperature::registerKeywords(Keywords& keys ){
+ ReweightBase::registerKeywords( keys );
+ keys.add("compulsory","REWEIGHT_TEMP","reweight data from a trajectory at one temperature and output the probability "
+ "distribution at a second temperature. This is not possible during postprocessing.");
+}
+
+ReweightTemperature::ReweightTemperature(const ActionOptions&ao):
+Action(ao),
+ReweightBase(ao)
+{
+ plumed.getAtoms().setCollectEnergy(true);
+ parse("REWEIGHT_TEMP",rtemp);
+ log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp);
+ rtemp*=plumed.getAtoms().getKBoltzmann();
+
+ retrieveAllBiases( "bias", biases );
+}
+
+double ReweightTemperature::getLogWeight() const {
+ // Retrieve the bias
+ double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
+ double energy=plumed.getAtoms().getEnergy()+bias;
+ return -( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias);
+}
+
+}
+}
diff --git a/src/core/ActionWithArguments.h b/src/core/ActionWithArguments.h
index 1321e6745..a881fd917 100644
--- a/src/core/ActionWithArguments.h
+++ b/src/core/ActionWithArguments.h
@@ -57,7 +57,7 @@ public:
/// Return a pointer to specific argument
Value* getPntrToArgument( const unsigned n );
/// Returns the number of arguments
- unsigned getNumberOfArguments() const ;
+ virtual unsigned getNumberOfArguments() const ;
/// Takes the difference taking into account pbc for arg i
double difference(int, double, double) const;
/// Takes one value and brings it back into the pbc of argument i
diff --git a/src/gridtools/ActionWithGrid.cpp b/src/gridtools/ActionWithGrid.cpp
index 464fd5c33..65ad56a62 100644
--- a/src/gridtools/ActionWithGrid.cpp
+++ b/src/gridtools/ActionWithGrid.cpp
@@ -20,116 +20,52 @@
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "ActionWithGrid.h"
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
namespace PLMD {
namespace gridtools {
void ActionWithGrid::registerKeywords( Keywords& keys ){
- Action::registerKeywords( keys );
- ActionPilot::registerKeywords( keys );
- ActionAtomistic::registerKeywords( keys );
- ActionWithArguments::registerKeywords( keys );
- ActionWithVessel::registerKeywords( keys );
- keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the grid");
- keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data. The default value "
- "of 0 implies that all the data will be used and that the grid will never be cleared");
+ vesselbase::ActionWithAveraging::registerKeywords( keys );
keys.add("compulsory","BANDWIDTH","the bandwidths for kernel density esimtation");
keys.add("compulsory","KERNEL","gaussian","the kernel function you are using. More details on the kernels available "
"in plumed plumed can be found in \\ref kernelfunctions.");
- keys.addFlag("UNORMALIZED",false,"output the unormalized density on the grid.");
}
ActionWithGrid::ActionWithGrid( const ActionOptions& ao):
Action(ao),
-ActionPilot(ao),
-ActionAtomistic(ao),
-ActionWithArguments(ao),
-ActionWithVessel(ao),
-requiresNorm(false),
+ActionWithAveraging(ao),
mygrid(NULL)
{
- if( keywords.exists("CLEAR") ){
- parse("CLEAR",clearstride);
- if( clearstride>0 ){
- if( clearstride%getStride()!=0 ) error("CLEAR parameter must be a multiple of STRIDE");
- log.printf(" clearing grid every %d steps \n",clearstride);
- }
- }
}
void ActionWithGrid::createGrid( const std::string& type, const std::string& inputstr ){
// Start creating the input for the grid
- std::string vstring = inputstr + " ";
+ std::string vstring = inputstr;
if( keywords.exists("KERNEL") ){
- vstring += getKeyword("KERNEL") + " " + getKeyword("BANDWIDTH");
+ std::string kstring; parse("KERNEL",kstring);
+ if( kstring=="DISCRETE" ) vstring += " KERNEL=" + kstring;
+ else vstring += " KERNEL=" + kstring + " " + getKeyword("BANDWIDTH");
}
- bool sumflag; parseFlag("UNORMALIZED",sumflag);
- if( sumflag ) vstring += " UNORMALIZED";
- if( clearstride>0 ) vstring += " NOMEMORY";
vesselbase::VesselOptions da("mygrid","",-1,vstring,this);
Keywords keys; gridtools::AverageOnGrid::registerKeywords( keys );
vesselbase::VesselOptions dar( da, keys );
if( type=="histogram" ){
- mygrid = new HistogramOnGrid(dar); requiresNorm=true; mygrid->setNorm(0);
+ mygrid = new HistogramOnGrid(dar);
} else if( type=="average" ){
- mygrid = new AverageOnGrid(dar); requiresNorm=false;
+ mygrid = new AverageOnGrid(dar);
} else if( type=="grid" ){
- mygrid = new GridVessel(dar); requiresNorm=false;
+ mygrid = new GridVessel(dar);
} else {
plumed_merror("no way to create grid of type " + type );
}
}
-void ActionWithGrid::finishGridSetup(){
- addVessel( mygrid );
-}
-
-void ActionWithGrid::lockRequests(){
- ActionAtomistic::lockRequests();
- ActionWithArguments::lockRequests();
-}
-
-void ActionWithGrid::unlockRequests(){
- ActionAtomistic::unlockRequests();
- ActionWithArguments::unlockRequests();
-}
-
-void ActionWithGrid::calculateNumericalDerivatives( PLMD::ActionWithValue* vv ){
- error("not possible to compute numerical derivatives for this action");
-}
-
-void ActionWithGrid::clearGrid(){
- mygrid->clear();
-}
-
-void ActionWithGrid::update(){
- if( getStep()==0 || !onStep() ) return;
- // Clear if it is time to reset
- if( mygrid ){
- if( mygrid->wasreset() ) clearGrid();
- }
- // Run all the tasks (if required
- if( prepareForTasks() ) runAllTasks();
- // This runs the jobs if we do not have a task list
- else performGridOperations( true );
- // Update the norm if this is necessary
- if( mygrid && requiresNorm ) mygrid->setNorm( 1 + mygrid->getNorm() );
- // Finish the tasks
- finishTaskSet();
- // By resetting here we are ensuring that the grid will be cleared at the start of the next step
- if( mygrid ){
- if( getStride()==0 || (clearstride>0 && getStep()%clearstride==0) ) mygrid->reset();
- }
-}
-
void ActionWithGrid::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
// Set the weight of this point
- myvals.setValue( 0, 1.0 ); compute( current, myvals );
-}
-
-void ActionWithGrid::performGridOperations( const bool& from_update ){
- plumed_error();
+ myvals.setValue( 0, cweight ); compute( current, myvals );
}
}
diff --git a/src/gridtools/ActionWithGrid.h b/src/gridtools/ActionWithGrid.h
index 604ad38ac..32d53e30d 100644
--- a/src/gridtools/ActionWithGrid.h
+++ b/src/gridtools/ActionWithGrid.h
@@ -22,53 +22,30 @@
#ifndef __PLUMED_gridtools_ActionWithGrid_h
#define __PLUMED_gridtools_ActionWithGrid_h
-#include "core/ActionPilot.h"
-#include "core/ActionAtomistic.h"
-#include "core/ActionWithArguments.h"
-#include "vesselbase/ActionWithVessel.h"
+#include "vesselbase/ActionWithAveraging.h"
#include "AverageOnGrid.h"
namespace PLMD {
namespace gridtools {
-class ActionWithGrid :
- public ActionPilot,
- public ActionAtomistic,
- public ActionWithArguments,
- public vesselbase::ActionWithVessel
-{
+class ActionWithGrid : public vesselbase::ActionWithAveraging {
private:
-/// Are we required to keep track of a normalization constant
- bool requiresNorm;
-/// The frequency with which to clear the grid
- unsigned clearstride;
/// The total number of bins
std::vector<unsigned> nbins;
/// The spacing between grid points
std::vector<double> gspacing;
-protected:
+/// The weights we are going to use for reweighting
+ std::vector<Value*> weights;
+protected:
/// The grid vessel
GridVessel* mygrid;
/// Read in stuff that is specifically for the grid and create it
void createGrid( const std::string& type, const std::string& inputstr );
-/// Finish the setup of the grid
- void finishGridSetup();
public:
static void registerKeywords( Keywords& keys );
explicit ActionWithGrid( const ActionOptions& );
- void lockRequests();
- void unlockRequests();
- void calculateNumericalDerivatives(PLMD::ActionWithValue*);
- void calculate(){}
- void apply(){}
- void update();
- virtual void clearGrid();
- virtual bool prepareForTasks() = 0;
- virtual void finishTaskSet() {};
void performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const ;
virtual void compute( const unsigned& current, MultiValue& myvals ) const = 0;
-/// This is used to perform all the tasks if we are not using the task list
- virtual void performGridOperations( const bool& from_update );
};
}
diff --git a/src/gridtools/ActionWithInputGrid.cpp b/src/gridtools/ActionWithInputGrid.cpp
index d57ee97ae..8d2569c39 100644
--- a/src/gridtools/ActionWithInputGrid.cpp
+++ b/src/gridtools/ActionWithInputGrid.cpp
@@ -58,22 +58,21 @@ ingrid(NULL)
log.printf(" using %uth component of grid calculated by action %s \n",mycomp,mves->getLabel().c_str() );
}
-bool ActionWithInputGrid::prepareForTasks(){
+void ActionWithInputGrid::clearAverage(){
+ mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), mygrid->getNbin(), mygrid->getGridSpacing() );
+ ActionWithAveraging::clearAverage();
+}
+
+void ActionWithInputGrid::prepareForAveraging(){
if( checkAllActive() ){
for(unsigned i=0;i<ingrid->getNumberOfPoints();++i){
if( ingrid->inactive(i) ) error("if FIND_CONTOUR is used with BUFFER option then other actions cannot be performed with grid");
}
}
- return (getFullNumberOfTasks()>0);
-}
-
-void ActionWithInputGrid::runFinalJobs(){
- if( getStride()>0 ) return;
- performGridOperations( false );
}
-void ActionWithInputGrid::performGridOperations( const bool& from_update ){
- plumed_assert( !from_update ); prepareForTasks(); runAllTasks();
+void ActionWithInputGrid::performOperations( const bool& from_update ){
+ prepareForAveraging(); runAllTasks();
}
}
diff --git a/src/gridtools/ActionWithInputGrid.h b/src/gridtools/ActionWithInputGrid.h
index 433aa4019..42823df59 100644
--- a/src/gridtools/ActionWithInputGrid.h
+++ b/src/gridtools/ActionWithInputGrid.h
@@ -40,10 +40,10 @@ protected:
public:
static void registerKeywords( Keywords& keys );
explicit ActionWithInputGrid(const ActionOptions&ao);
- bool prepareForTasks();
- void runFinalJobs();
+ virtual void clearAverage();
+ virtual void prepareForAveraging();
virtual bool checkAllActive() const { return true; }
- virtual void performGridOperations( const bool& from_update );
+ virtual void performOperations( const bool& from_update );
};
inline
diff --git a/src/gridtools/AverageOnGrid.cpp b/src/gridtools/AverageOnGrid.cpp
index 3167b53b7..0f8933f1c 100644
--- a/src/gridtools/AverageOnGrid.cpp
+++ b/src/gridtools/AverageOnGrid.cpp
@@ -49,23 +49,20 @@ void AverageOnGrid::accumulate( const unsigned& ipoint, const double& weight, co
}
double AverageOnGrid::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
- if( unormalised ) return data[nper*ipoint + jelement];
+ if( noAverage() ) return getDataElement( nper*ipoint + jelement);
- if( noderiv ) return data[nper*ipoint+jelement] / data[nper*(1+ipoint) - 1];
+ if( jelement>=(nper-(dimension+1)) ) return getDataElement( nper*ipoint + jelement );
+
+ if( noderiv ) return getDataElement( nper*ipoint+jelement ) / getDataElement( nper*(1+ipoint) - 1);
double rdenom = 1.0;
- if( fabs(data[nper*(ipoint+1) -(dimension+1)])>epsilon ) rdenom = 1. / data[nper*(ipoint+1) - (dimension+1)];
+ if( fabs(getDataElement( nper*(ipoint+1) -(dimension+1) ))>epsilon ) rdenom = 1. / getDataElement( nper*(ipoint+1) - (dimension+1) );
unsigned jderiv = jelement%(1+dimension);
- if( jderiv==0 ) return rdenom*data[nper*ipoint+jelement];
+ if( jderiv==0 ) return rdenom*getDataElement( nper*ipoint+jelement );
unsigned jfloor = std::floor( jelement / (1+dimension) );
- return rdenom*data[nper*ipoint+jelement] - rdenom*rdenom*data[nper*ipoint+jfloor]*data[nper*(ipoint+1) - (dimension+1) + jderiv ];
-}
-
-double AverageOnGrid::getGridElementForPrint( const unsigned& ipoint, const unsigned& jelement ) const {
- plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
- return data[nper*ipoint + jelement];
+ return rdenom*getDataElement( nper*ipoint+jelement ) - rdenom*rdenom*getDataElement(nper*ipoint+jfloor)*getDataElement(nper*(ipoint+1) - (dimension+1) + jderiv);
}
}
diff --git a/src/gridtools/AverageOnGrid.h b/src/gridtools/AverageOnGrid.h
index 5229e0e52..7288526fa 100644
--- a/src/gridtools/AverageOnGrid.h
+++ b/src/gridtools/AverageOnGrid.h
@@ -33,7 +33,6 @@ public:
explicit AverageOnGrid( const vesselbase::VesselOptions& da );
void accumulate( const unsigned& ipoint, const double& weight, const double& dens, const std::vector<double>& der, std::vector<double>& buffer ) const ;
double getGridElement( const unsigned& ipoint, const unsigned& jelement ) const ;
- double getGridElementForPrint( const unsigned& ipoint, const unsigned& jelement ) const ;
unsigned getNumberOfComponents() const ;
};
diff --git a/src/gridtools/ContourFindingBase.cpp b/src/gridtools/ContourFindingBase.cpp
index be63da700..7d3d2ccad 100644
--- a/src/gridtools/ContourFindingBase.cpp
+++ b/src/gridtools/ContourFindingBase.cpp
@@ -80,7 +80,7 @@ mydata(NULL)
}
}
-void ContourFindingBase::finishTaskSet(){
+void ContourFindingBase::finishAveraging(){
// And this looks after the output
if( mydata ){
std::vector<double> point( 1 + ingrid->getDimension() );
diff --git a/src/gridtools/ContourFindingBase.h b/src/gridtools/ContourFindingBase.h
index 90a6c0d1e..2dbf0dbbc 100644
--- a/src/gridtools/ContourFindingBase.h
+++ b/src/gridtools/ContourFindingBase.h
@@ -55,7 +55,7 @@ public:
/// Number of quantities is the number of points in each point on the grid
virtual unsigned getNumberOfQuantities() const { return 1 + ingrid->getDimension(); }
/// This does output if needs be
- void finishTaskSet();
+ void finishAveraging();
};
inline
diff --git a/src/gridtools/ConvertToFES.cpp b/src/gridtools/ConvertToFES.cpp
index d91e656eb..1712f159c 100644
--- a/src/gridtools/ConvertToFES.cpp
+++ b/src/gridtools/ConvertToFES.cpp
@@ -54,7 +54,8 @@ PLUMED_REGISTER_ACTION(ConvertToFES,"CONVERT_TO_FES")
void ConvertToFES::registerKeywords( Keywords& keys ){
ActionWithInputGrid::registerKeywords( keys );
keys.add("optional","TEMP","the temperature at which you are operating");
- keys.remove("STRIDE"); keys.remove("KERNEL"); keys.remove("BANDWIDTH"); keys.remove("CLEAR");
+ keys.remove("STRIDE"); keys.remove("KERNEL"); keys.remove("BANDWIDTH");
+ keys.remove("LOGWEIGHTS"); keys.remove("CLEAR");
}
ConvertToFES::ConvertToFES(const ActionOptions&ao):
@@ -68,7 +69,7 @@ ActionWithInputGrid(ao)
if( ingrid->noDerivatives() ) mygrid->setNoDerivatives();
std::vector<double> fspacing;
mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), ingrid->getNbin(), fspacing);
- finishGridSetup(); resizeFunctions();
+ setAveragingAction( mygrid, true );
simtemp=0.; parse("TEMP",simtemp);
if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
diff --git a/src/gridtools/DumpCube.cpp b/src/gridtools/DumpCube.cpp
index 871bfc68a..67109b181 100644
--- a/src/gridtools/DumpCube.cpp
+++ b/src/gridtools/DumpCube.cpp
@@ -20,7 +20,7 @@
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
#include "core/ActionRegister.h"
-#include "ActionWithInputGrid.h"
+#include "GridPrintingBase.h"
#include "tools/OFile.h"
//+PLUMEDOC GRIDANALYSIS DUMPCUBE
@@ -35,46 +35,42 @@ Output a three dimensional grid using the Gaussian cube file format.
namespace PLMD {
namespace gridtools {
-class DumpCube : public ActionWithInputGrid {
+class DumpCube : public GridPrintingBase {
private:
- std::string fmt, filename;
+ unsigned mycomp;
public:
static void registerKeywords( Keywords& keys );
explicit DumpCube(const ActionOptions&ao);
- void performGridOperations( const bool& from_update );
- unsigned getNumberOfDerivatives(){ return 0; }
- void compute( const unsigned& , MultiValue& ) const {}
- bool isPeriodic(){ return false; }
+ void printGrid( OFile& ofile ) const ;
};
PLUMED_REGISTER_ACTION(DumpCube,"DUMPCUBE")
void DumpCube::registerKeywords( Keywords& keys ){
- ActionWithInputGrid::registerKeywords( keys );
- keys.add("compulsory","FILE","density","the file on which to write the grid.");
- keys.add("optional","FMT","the format that should be used to output real numbers");
+ GridPrintingBase::registerKeywords( keys );
+ keys.add("optional","COMPONENT","if your input is a vector field use this to specifiy the component of the input vector field for which you wish to output");
}
DumpCube::DumpCube(const ActionOptions&ao):
Action(ao),
-ActionWithInputGrid(ao),
-fmt("%f")
+GridPrintingBase(ao)
{
+ fmt = fmt + " ";
if( ingrid->getDimension()!=3 ) error("cannot print cube file if grid does not contain three dimensional data");
- parse("FMT",fmt); fmt=fmt +" "; parse("FILE",filename);
- if(filename.length()==0) error("name out output file was not specified");
- log.printf(" outputting grid to file named %s with format %s \n",filename.c_str(), fmt.c_str() );
+ if( ingrid->getNumberOfComponents()==1 ){
+ mycomp=0;
+ } else {
+ int tcomp=-1; parse("COMPONENT",tcomp);
+ if( tcomp<0 ) error("component of vector field was not specified - use COMPONENT keyword");
+ mycomp=tcomp*(1+ingrid->getDimension()); if( ingrid->noDerivatives() ) mycomp=tcomp;
+ log.printf(" using %dth component of grid \n",tcomp );
+ }
checkRead();
}
-void DumpCube::performGridOperations( const bool& from_update ){
- if( !from_update && getStride()>0 ) return ;
- OFile ofile; ofile.link(*this);
- if( from_update ) ofile.setBackupString("analysis");
- ofile.open( filename );
-
+void DumpCube::printGrid( OFile& ofile ) const {
double lunit = ingrid->getCubeUnits();
ofile.printf("PLUMED CUBE FILE\n");
@@ -90,14 +86,12 @@ void DumpCube::performGridOperations( const bool& from_update ){
for(pp[0]=0;pp[0]<nbin[0];++pp[0]){
for(pp[1]=0;pp[1]<nbin[1];++pp[1]){
for(pp[2]=0;pp[2]<nbin[2];++pp[2]){
- ofile.printf(fmt.c_str(), getFunctionValue(pp) );
+ ofile.printf(fmt.c_str(), ingrid->getGridElement( ingrid->getIndex(pp), mycomp ) );
if(pp[2]%6==5) ofile.printf("\n");
}
ofile.printf("\n");
}
}
-
- ofile.close();
}
}
diff --git a/src/gridtools/DumpGrid.cpp b/src/gridtools/DumpGrid.cpp
index c2ebec60a..150378ca1 100644
--- a/src/gridtools/DumpGrid.cpp
+++ b/src/gridtools/DumpGrid.cpp
@@ -19,10 +19,8 @@
You should have received a copy of the GNU Lesser General Public License
along with plumed. If not, see <http://www.gnu.org/licenses/>.
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
-#include "ActionWithInputGrid.h"
+#include "GridPrintingBase.h"
#include "core/ActionRegister.h"
-#include "AverageOnGrid.h"
-#include "tools/IFile.h"
#include "tools/OFile.h"
namespace PLMD {
@@ -37,50 +35,27 @@ Output the function on the grid to a file with the PLUMED grid format.
*/
//+ENDPLUMEDOC
-class DumpGrid : public ActionWithInputGrid {
-private:
- bool printav;
- std::string fmt, filename;
+class DumpGrid : public GridPrintingBase {
public:
static void registerKeywords( Keywords& keys );
explicit DumpGrid(const ActionOptions&ao);
- void performGridOperations( const bool& from_update );
- unsigned getNumberOfDerivatives(){ return 0; }
- void compute( const unsigned& , MultiValue& ) const {}
- bool isPeriodic(){ return false; }
+ void printGrid( OFile& ofile ) const ;
};
PLUMED_REGISTER_ACTION(DumpGrid,"DUMPGRID")
void DumpGrid::registerKeywords( Keywords& keys ){
- ActionWithInputGrid::registerKeywords( keys );
- keys.add("compulsory","FILE","density","the file on which to write the grid.");
- keys.add("optional","FMT","the format that should be used to output real numbers");
- keys.addFlag("PRINT_AVERAGE",false,"if your input is a \\ref MULTICOLVARDENS output the average in the grid file.");
+ GridPrintingBase::registerKeywords( keys );
}
DumpGrid::DumpGrid(const ActionOptions&ao):
Action(ao),
-ActionWithInputGrid(ao),
-printav(false),
-fmt("%f")
+GridPrintingBase(ao)
{
- AverageOnGrid* myav = dynamic_cast<AverageOnGrid*>( ingrid );
- if( myav ) parseFlag("PRINT_AVERAGE",printav);
-
- parse("FMT",fmt); fmt=" "+fmt; parse("FILE",filename);
- if(filename.length()==0) error("name out output file was not specified");
- log.printf(" outputting grid to file named %s with format %s \n",filename.c_str(), fmt.c_str() );
- checkRead();
+ fmt = " " + fmt; checkRead();
}
-void DumpGrid::performGridOperations( const bool& from_update ){
- if( !from_update && getStride()>0 ) return ;
-
- OFile ofile; ofile.link(*this);
- if ( from_update ) ofile.setBackupString("analysis");
- ofile.open( filename );
-
+void DumpGrid::printGrid( OFile& ofile ) const {
ofile.addConstantField("normalisation");
for(unsigned i=0;i<ingrid->getDimension();++i){
ofile.addConstantField("min_" + ingrid->getComponentName(i) );
@@ -94,7 +69,7 @@ void DumpGrid::performGridOperations( const bool& from_update ){
for(unsigned i=0;i<ingrid->getNumberOfPoints();++i){
ingrid->getIndices( i, ind );
if(i>0 && ingrid->getDimension()==2 && ind[ingrid->getDimension()-2]==0) ofile.printf("\n");
- ofile.fmtField(fmt); ofile.printField("normalisation", ingrid->norm );
+ ofile.fmtField(fmt); ofile.printField("normalisation", ingrid->getNorm() );
for(unsigned j=0;j<ingrid->getDimension();++j){
ofile.printField("min_" + ingrid->getComponentName(j), ingrid->getMin()[j] );
ofile.printField("max_" + ingrid->getComponentName(j), ingrid->getMax()[j] );
@@ -105,19 +80,11 @@ void DumpGrid::performGridOperations( const bool& from_update ){
// Retrieve and print the grid coordinates
ingrid->getGridPointCoordinates(i, xx );
for(unsigned j=0;j<ingrid->getDimension();++j){ ofile.fmtField(fmt); ofile.printField(ingrid->getComponentName(j),xx[j]); }
- if( printav ){
- unsigned nnorm=ingrid->dimension+1; if( ingrid->noderiv ) nnorm=1;
- for(unsigned j=0;j<ingrid->getNumberOfQuantities()-nnorm;++j){
- ofile.fmtField(fmt); ofile.printField(ingrid->arg_names[ingrid->dimension+j], ingrid->getGridElement( i, j ) );
- }
- } else {
- for(unsigned j=0;j<ingrid->getNumberOfQuantities();++j){
- ofile.fmtField(fmt); ofile.printField(ingrid->arg_names[ingrid->dimension+j], ingrid->getGridElementForPrint( i, j ) );
- }
+ for(unsigned j=0;j<ingrid->getNumberOfQuantities();++j){
+ ofile.fmtField(fmt); ofile.printField(ingrid->arg_names[ingrid->dimension+j], ingrid->getGridElement( i, j ) );
}
ofile.printField();
}
- ofile.close();
}
}
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index bb7058a85..a60c2ac59 100644
--- a/src/gridtools/FindContour.cpp
+++ b/src/gridtools/FindContour.cpp
@@ -42,10 +42,10 @@ public:
static void registerKeywords( Keywords& keys );
explicit FindContour(const ActionOptions&ao);
bool checkAllActive() const { return gbuffer==0; }
- bool prepareForTasks();
+ void prepareForAveraging();
bool isPeriodic(){ return false; }
void compute( const unsigned& current, MultiValue& myvals ) const ;
- void finishTaskSet();
+ void finishAveraging();
};
PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
@@ -67,7 +67,7 @@ firsttime(true)
checkRead();
}
-bool FindContour::prepareForTasks(){
+void FindContour::prepareForAveraging(){
// Create a task list if first time
if( firsttime ){
for(unsigned i=0;i<ingrid->getDimension()*ingrid->getNumberOfPoints();++i) addTaskToList( i );
@@ -106,7 +106,7 @@ bool FindContour::prepareForTasks(){
else ind[j]-=1;
}
}
- lockContributors(); return true;
+ lockContributors();
}
void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
@@ -126,8 +126,8 @@ void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
for(unsigned i=0;i<ingrid->getDimension();++i) myvals.setValue( 1+i, point[i] );
}
-void FindContour::finishTaskSet(){
- ContourFindingBase::finishTaskSet();
+void FindContour::finishAveraging(){
+ ContourFindingBase::finishAveraging();
// And update the list of active grid points
if( gbuffer>0 ){
std::vector<unsigned> neighbours; unsigned num_neighbours;
diff --git a/src/gridtools/FindContourSurface.cpp b/src/gridtools/FindContourSurface.cpp
index dce11fed5..cc10f077a 100644
--- a/src/gridtools/FindContourSurface.cpp
+++ b/src/gridtools/FindContourSurface.cpp
@@ -47,9 +47,10 @@ public:
explicit FindContourSurface(const ActionOptions&ao);
unsigned getNumberOfQuantities() const { return 2; }
bool checkAllActive() const { return gbuffer==0; }
- bool prepareForTasks();
+ void clearAverage();
+ void prepareForAveraging();
void compute( const unsigned& current, MultiValue& myvals ) const ;
- void finishTaskSet();
+ void finishAveraging();
};
PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
@@ -94,18 +95,23 @@ ones(ingrid->getDimension(),1)
for(unsigned i=1;i<gdirs.size();++i){
if( ingrid->isPeriodic(gdirs[i]) ) vstring+=",T"; else vstring+=",F";
}
- createGrid( "grid", vstring ); mygrid->setNoDerivatives(); finishGridSetup();
+ createGrid( "grid", vstring ); mygrid->setNoDerivatives();
+ setAveragingAction( mygrid, true );
}
-bool FindContourSurface::prepareForTasks(){
+void FindContourSurface::clearAverage(){
// Set the boundaries of the output grid
std::vector<double> fspacing; std::vector<unsigned> snbins( ingrid->getDimension()-1 );
std::vector<std::string> smin( ingrid->getDimension()-1 ), smax( ingrid->getDimension()-1 );
for(unsigned i=0;i<gdirs.size();++i){
smin[i]=ingrid->getMin()[gdirs[i]]; smax[i]=ingrid->getMax()[gdirs[i]];
snbins[i]=ingrid->getNbin()[gdirs[i]];
- }
+ }
mygrid->setBounds( smin, smax, snbins, fspacing); resizeFunctions();
+ ActionWithAveraging::clearAverage();
+}
+
+void FindContourSurface::prepareForAveraging(){
// Create a task list if first time
if( firsttime ){
std::vector<unsigned> find( ingrid->getDimension() );
@@ -124,11 +130,11 @@ bool FindContourSurface::prepareForTasks(){
direction.resize( ingrid->getDimension(), 0 );
direction[dir_n] = 0.999999999*ingrid->getGridSpacing()[dir_n];
}
- firsttime=false; return true;
+ firsttime=false;
}
-void FindContourSurface::finishTaskSet(){
- ContourFindingBase::finishTaskSet();
+void FindContourSurface::finishAveraging(){
+ ContourFindingBase::finishAveraging();
// And update the list of active grid points
if( gbuffer>0 ){
std::vector<double> dx( ingrid->getGridSpacing() );
diff --git a/src/gridtools/FourierTransform.cpp b/src/gridtools/FourierTransform.cpp
index a6ec20e2b..9a363c9de 100644
--- a/src/gridtools/FourierTransform.cpp
+++ b/src/gridtools/FourierTransform.cpp
@@ -76,9 +76,9 @@ public:
static void registerKeywords( Keywords& keys );
explicit FourierTransform(const ActionOptions&ao);
#ifndef __PLUMED_HAS_FFTW
- void performOperationsWithGrid( const bool& from_update ){}
+ void performOperations( const bool& from_update ){}
#else
- void performOperationsWithGrid( const bool& from_update );
+ void performOperations( const bool& from_update );
#endif
unsigned getNumberOfDerivatives(){ return 0; }
void compute( const unsigned& , MultiValue& ) const {}
@@ -173,7 +173,7 @@ fourier_params(2)
}
#ifdef __PLUMED_HAS_FFTW
-void FourierTransform::performOperationsWithGrid( const bool& from_update ){
+void FourierTransform::performOperations( const bool& from_update ){
// Spacing of the real grid
std::vector<double> g_spacing ( ingrid->getGridSpacing() );
diff --git a/src/gridtools/GridPrintingBase.cpp b/src/gridtools/GridPrintingBase.cpp
new file mode 100644
index 000000000..c2d253c04
--- /dev/null
+++ b/src/gridtools/GridPrintingBase.cpp
@@ -0,0 +1,76 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2014,2015 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "GridPrintingBase.h"
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
+#include "vesselbase/ActionWithVessel.h"
+
+namespace PLMD {
+namespace gridtools {
+
+void GridPrintingBase::registerKeywords( Keywords& keys ){
+ Action::registerKeywords( keys ); ActionPilot::registerKeywords( keys );
+ keys.add("compulsory","GRID","the action that creates the grid you would like to output");
+ keys.add("compulsory","STRIDE","0","the frequency with which the grid should be output to the file. The default "
+ "value of 0 ensures that the grid is only output at the end of the trajectory");
+ keys.add("compulsory","FILE","density","the file on which to write the grid.");
+ keys.add("optional","FMT","the format that should be used to output real numbers");
+}
+
+GridPrintingBase::GridPrintingBase(const ActionOptions&ao):
+Action(ao),
+ActionPilot(ao),
+fmt("%f")
+{
+ std::string mlab; parse("GRID",mlab);
+ vesselbase::ActionWithVessel* mves= plumed.getActionSet().selectWithLabel<vesselbase::ActionWithVessel*>(mlab);
+ if(!mves) error("action labelled " + mlab + " does not exist or does not have vessels");
+ addDependency(mves);
+
+ for(unsigned i=0;i<mves->getNumberOfVessels();++i){
+ ingrid=dynamic_cast<GridVessel*>( mves->getPntrToVessel(i) );
+ if( ingrid ) break;
+ }
+ if( !ingrid ) error("input action does not calculate a grid");
+
+ parse("FMT",fmt); parse("FILE",filename);
+ if(filename.length()==0) error("name out output file was not specified");
+ log.printf(" outputting grid calculated by action %s to file named %s with format %s \n",mves->getLabel().c_str(),filename.c_str(), fmt.c_str() );
+}
+
+void GridPrintingBase::update(){
+ if( getStep()==0 || getStride()==0 ) return ;
+
+ OFile ofile; ofile.link(*this);
+ ofile.setBackupString("analysis");
+ ofile.open( filename ); printGrid( ofile ); ofile.close();
+}
+
+void GridPrintingBase::runFinalJobs(){
+ if( getStride()>0 ) return;
+
+ OFile ofile; ofile.link(*this);
+ ofile.open( filename ); printGrid( ofile ); ofile.close();
+}
+
+}
+}
diff --git a/src/gridtools/GridPrintingBase.h b/src/gridtools/GridPrintingBase.h
new file mode 100644
index 000000000..b6444b022
--- /dev/null
+++ b/src/gridtools/GridPrintingBase.h
@@ -0,0 +1,48 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2014,2015 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_gridtools_GridPrintingBase_h
+#define __PLUMED_gridtools_GridPrintingBase_h
+
+#include "core/ActionPilot.h"
+#include "GridVessel.h"
+#include "tools/OFile.h"
+
+namespace PLMD {
+namespace gridtools {
+
+class GridPrintingBase : public ActionPilot {
+protected:
+ GridVessel* ingrid;
+ std::string fmt, filename;
+public:
+ static void registerKeywords( Keywords& keys );
+ explicit GridPrintingBase(const ActionOptions&ao);
+ void calculate(){}
+ void apply(){}
+ void update();
+ void runFinalJobs();
+ virtual void printGrid( OFile& ofile ) const=0;
+};
+
+}
+}
+#endif
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 2f7288590..dfd383189 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -27,20 +27,16 @@ namespace PLMD {
namespace gridtools {
void GridVessel::registerKeywords( Keywords& keys ){
- Vessel::registerKeywords( keys );
+ AveragingVessel::registerKeywords( keys );
keys.add("compulsory","COMPONENTS","the names of the components in the vector");
keys.add("compulsory","COORDINATES","the names of the coordinates of the grid");
keys.add("compulsory","PBC","is the grid periodic in each direction or not");
- keys.addFlag("NOMEMORY",false,"should the data in the grid be deleted after print/analysis");
- keys.addFlag("UNORMALIZED",false,"don't normalise grid");
}
GridVessel::GridVessel( const vesselbase::VesselOptions& da ):
-Vessel(da),
+AveragingVessel(da),
bounds_set(false),
cube_units(1.0),
-wascleared(true),
-norm(0.0),
noderiv(false)
{
std::vector<std::string> compnames; parseVector("COMPONENTS",compnames);
@@ -48,8 +44,6 @@ noderiv(false)
dimension=coordnames.size();
std::vector<std::string> spbc( dimension ); parseVector("PBC",spbc);
str_min.resize( dimension); str_max.resize( dimension );
- parseFlag("NOMEMORY",nomemory); parseFlag("UNORMALIZED",unormalised);
- foundprint=nomemory;
unsigned n=0; nper=compnames.size()*( 1 + coordnames.size() );
arg_names.resize( coordnames.size() + compnames.size()*( 1 + coordnames.size() ) );
@@ -90,9 +84,10 @@ void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vec
Tools::convert( str_min[i], min[i] );
Tools::convert( str_max[i], max[i] );
if( spacing.size()==dimension && binsin.size()==dimension ){
- unsigned spc = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
- if( spc>binsin[i] ) nbin[i]=spc;
- else nbin[i]=binsin[i];
+ double range = max[i] - min[i]; unsigned spc = std::floor( range / spacing[i]);
+ // This check ensures that nbins is set correctly if spacing is set the same as the number of bins
+ if( fabs( binsin[i]*spacing[i]-range )>epsilon ) spc += 1;
+ if( spc>binsin[i] ) nbin[i]=spc; else nbin[i]=binsin[i];
} else if( binsin.size()==dimension ) nbin[i]=binsin[i];
else if( spacing.size()==dimension ) nbin[i] = std::floor(( max[i] - min[i] ) / spacing[i]) + 1;
else plumed_error();
@@ -101,6 +96,7 @@ void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vec
stride[i]=npoints;
npoints*=nbin[i];
}
+ resize(); // Always resize after setting new bounds as grid size may have have changed
}
std::string GridVessel::description(){
@@ -123,11 +119,8 @@ std::string GridVessel::description(){
void GridVessel::resize(){
plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
- resizeBuffer( getNumberOfBufferPoints()*nper );
- if( data.size()!=npoints*nper ){
- data.resize( npoints*nper, 0 );
- active.resize( npoints, true );
- }
+ resizeBuffer( getNumberOfBufferPoints()*nper ); setDataSize( npoints*nper );
+ if( active.size()!=npoints) active.resize( npoints, true );
}
unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
@@ -198,16 +191,12 @@ void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned
double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
- return data[ nper*ipoint + jelement ];
-}
-
-double GridVessel::getGridElementForPrint( const unsigned& ipoint, const unsigned& jelement ) const {
- return getGridElement( ipoint, jelement );
+ return getDataElement( nper*ipoint + jelement );
}
void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ){
plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
- wascleared=false; data[ nper*ipoint + jelement ] = value;
+ setDataElement( nper*ipoint + jelement, value );
}
void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
@@ -215,10 +204,6 @@ void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::ve
for(unsigned i=0;i<nper;++i) buffer[bufstart + nper*current + i] += myvals.get(i+1);
}
-void GridVessel::finish( const std::vector<double>& buffer ){
- wascleared=false; for(unsigned i=0;i<data.size();++i) data[i]+=buffer[bufstart + i];
-}
-
double GridVessel::getGridElement( const std::vector<unsigned>& indices, const unsigned& jelement ) const {
return getGridElement( getIndex( indices ), jelement );
}
@@ -285,15 +270,6 @@ void GridVessel::getNeighbors( const std::vector<unsigned>& indices, const std::
}
}
-void GridVessel::reset(){
- wascleared=true;
-}
-
-void GridVessel::clear(){
- plumed_assert( wasreset() );
- norm=0; data.assign( data.size(), 0.0 );
-}
-
void GridVessel::setCubeUnits( const double& units ){
cube_units=units;
}
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 0cb5dfb09..83a215633 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -25,15 +25,13 @@
#include <string>
#include <cstring>
#include <vector>
-#include "vesselbase/Vessel.h"
+#include "vesselbase/AveragingVessel.h"
namespace PLMD {
namespace gridtools {
-class GridVessel : public vesselbase::Vessel {
+class GridVessel : public vesselbase::AveragingVessel {
friend class ActionWithInputGrid;
-friend class AverageOnGrid;
-friend class GridFunction;
friend class DumpGrid;
private:
/// Have the minimum and maximum for the grid been set
@@ -53,18 +51,10 @@ private:
/// The grid point that was requested last by getGridPointCoordinates
unsigned currentGridPoint;
protected:
-/// The grid was recently cleared and bounds can be set
- bool wascleared;
/// Do we have derivatives
bool noderiv;
/// The names of the various columns in the grid file
std::vector<std::string> arg_names;
-/// The normalisation constant to use
- double norm;
-/// Are we deleting the data after print
- bool nomemory;
-/// Are we outputting unormalised data
- bool unormalised;
/// The number of pieces of information we are storing for each
/// point in the grid
unsigned nper;
@@ -78,8 +68,6 @@ protected:
unsigned dimension;
/// Which grid points are we actively accumulating
std::vector<bool> active;
-/// The grid with all the data values on it
- std::vector<double> data;
/// Convert a point in space the the correspoinding grid point
unsigned getIndex( const std::vector<double>& p ) const ;
public:
@@ -138,7 +126,6 @@ public:
double getCellVolume() const ;
/// Get the value of the ith grid element
virtual double getGridElement( const unsigned&, const unsigned& ) const ;
- virtual double getGridElementForPrint( const unsigned&, const unsigned& ) const ;
/// Get the set of points neighouring a particular location in space
void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
@@ -153,33 +140,20 @@ public:
double getGridExtent( const unsigned& i ) const ;
/// Copy data from the action into the grid
virtual void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
-/// Copy data from an accumulated buffer into the grid
- virtual void finish( const std::vector<double>& );
-/// Clear all the data stored on the grid
- virtual void clear();
-/// Reset the grid so that it is cleared at start of next time it is calculated
- virtual void reset();
/// This ensures that Gaussian cube fies are in correct units
void setCubeUnits( const double& units );
/// This ensures that Gaussian cube files are in correct units
double getCubeUnits() const ;
- virtual void switchOffNormalisation(){}
/// Return a string containing the input to the grid so we can clone it
std::string getInputString() const ;
/// Does this have derivatives
bool noDerivatives() const ;
/// Get the value and derivatives at a particular location using spline interpolation
double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ;
-/// Was the grid cleared on the last step
- bool wasreset() const ;
/// Deactivate all the grid points
void activateThesePoints( const std::vector<bool>& to_activate );
/// Is this point active
bool inactive( const unsigned& ip ) const ;
-/// Functions for dealing with normalisation constant
- void setNorm( const double& snorm );
- double getNorm() const ;
- bool applyForce( std::vector<double>& forces ){ return false; }
};
inline
@@ -234,27 +208,12 @@ bool GridVessel::noDerivatives() const {
return noderiv;
}
-inline
-bool GridVessel::wasreset() const {
- return wascleared;
-}
-
inline
bool GridVessel::inactive( const unsigned& ip ) const {
plumed_dbg_assert( ip<npoints );
return !active[ip];
}
-inline
-void GridVessel::setNorm( const double& snorm ){
- norm=snorm;
-}
-
-inline
-double GridVessel::getNorm() const {
- return norm;
-}
-
inline
const std::vector<unsigned>& GridVessel::getStride() const {
return stride;
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 36725555b..a6468de4b 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -46,6 +46,10 @@ discrete(false)
}
}
+bool HistogramOnGrid::noDiscreteKernels() const {
+ return !discrete;
+}
+
void HistogramOnGrid::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
const std::vector<unsigned>& nbins, const std::vector<double>& spacing ){
GridVessel::setBounds( smin, smax, nbins, spacing );
@@ -103,7 +107,6 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
} else {
std::vector<Value*> vv( getVectorOfValues() );
- unsigned ntot=nper*getNumberOfPoints();
double newval; std::vector<double> xx( dimension );
for(unsigned i=0;i<num_neigh;++i){
unsigned ineigh=neighbors[i];
@@ -126,19 +129,13 @@ void HistogramOnGrid::accumulate( const unsigned& ipoint, const double& weight,
void HistogramOnGrid::finish( const std::vector<double>& buffer ){
if( addOneKernelAtATime ){
- wascleared=false;
for(unsigned i=0;i<getAction()->getCurrentNumberOfActiveTasks();++i){
- for(unsigned j=0;j<nper;++j) data[nper*getAction()->getActiveTask(i)+j]+=buffer[bufstart+i*nper+j];
+ for(unsigned j=0;j<nper;++j) addDataElement( nper*getAction()->getActiveTask(i)+j, buffer[bufstart+i*nper+j] );
}
} else {
GridVessel::finish( buffer );
}
}
-double HistogramOnGrid::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
- if( unormalised ) return GridVessel::getGridElement( ipoint, jelement );
- return GridVessel::getGridElement( ipoint, jelement ) / norm;
-}
-
}
}
diff --git a/src/gridtools/HistogramOnGrid.h b/src/gridtools/HistogramOnGrid.h
index 5670118df..f26d07c1f 100644
--- a/src/gridtools/HistogramOnGrid.h
+++ b/src/gridtools/HistogramOnGrid.h
@@ -47,11 +47,11 @@ public:
void calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
void finish( const std::vector<double>& buffer );
virtual void accumulate( const unsigned& ipoint, const double& weight, const double& dens, const std::vector<double>& der, std::vector<double>& buffer ) const ;
- virtual double getGridElement( const unsigned& ipoint, const unsigned& jelement ) const ;
unsigned getNumberOfBufferPoints() const ;
KernelFunctions* getKernelAndNeighbors( std::vector<double>& point, unsigned& num_neigh, std::vector<unsigned>& neighbors ) const;
std::vector<Value*> getVectorOfValues() const ;
void addOneKernelEachTimeOnly(){ addOneKernelAtATime=true; }
+ bool noDiscreteKernels() const ;
};
inline
diff --git a/src/gridtools/InterpolateGrid.cpp b/src/gridtools/InterpolateGrid.cpp
index 16c8fdc58..763799c0c 100644
--- a/src/gridtools/InterpolateGrid.cpp
+++ b/src/gridtools/InterpolateGrid.cpp
@@ -37,18 +37,13 @@ namespace PLMD {
namespace gridtools {
class InterpolateGrid : public ActionWithInputGrid {
-private:
- std::vector<unsigned> nbin;
- std::vector<double> gspacing;
public:
static void registerKeywords( Keywords& keys );
explicit InterpolateGrid(const ActionOptions&ao);
unsigned getNumberOfDerivatives(){ return 0; }
unsigned getNumberOfQuantities() const ;
- void clearGrid();
void compute( const unsigned& current, MultiValue& myvals ) const ;
bool isPeriodic(){ return false; }
- bool onStep() const { return true; }
};
PLUMED_REGISTER_ACTION(InterpolateGrid,"INTERPOLATE_GRID")
@@ -57,7 +52,7 @@ void InterpolateGrid::registerKeywords( Keywords& keys ){
ActionWithInputGrid::registerKeywords( keys );
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.remove("STRIDE"); keys.remove("KERNEL"); keys.remove("BANDWIDTH"); keys.remove("CLEAR");
+ keys.remove("KERNEL"); keys.remove("BANDWIDTH");
}
InterpolateGrid::InterpolateGrid(const ActionOptions&ao):
@@ -68,16 +63,16 @@ ActionWithInputGrid(ao)
if( ingrid->noDerivatives() ) error("cannot interpolate a grid that does not have derivatives");
// Create the input from the old string
createGrid( "grid", "COMPONENTS=" + getLabel() + " " + ingrid->getInputString() );
- finishGridSetup();
- parseVector("GRID_BIN",nbin); parseVector("GRID_SPACING",gspacing);
+ std::vector<unsigned> nbin; parseVector("GRID_BIN",nbin);
+ std::vector<double> gspacing; parseVector("GRID_SPACING",gspacing);
if( nbin.size()!=ingrid->getDimension() && gspacing.size()!=ingrid->getDimension() ){
error("GRID_BIN or GRID_SPACING must be set");
}
// Need this for creation of tasks
mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), nbin, gspacing );
- resizeFunctions();
+ setAveragingAction( mygrid, true );
// Now create task list
for(unsigned i=0;i<mygrid->getNumberOfPoints();++i) addTaskToList(i);
@@ -91,10 +86,6 @@ unsigned InterpolateGrid::getNumberOfQuantities() const {
return 2 + ingrid->getDimension();
}
-void InterpolateGrid::clearGrid(){
- mygrid->setBounds( ingrid->getMin(), ingrid->getMax(), nbin, gspacing ); mygrid->clear();
-}
-
void InterpolateGrid::compute( const unsigned& current, MultiValue& myvals ) const {
std::vector<double> pos( mygrid->getDimension() ); mygrid->getGridPointCoordinates( current, pos );
std::vector<double> der( mygrid->getDimension() ); double val = getFunctionValueAndDerivatives( pos, der );
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index e7eece5c2..eae18b993 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/src/multicolvar/MultiColvarDensity.cpp
@@ -101,8 +101,8 @@ public:
unsigned getNumberOfQuantities() const ;
bool isPeriodic(){ return false; }
unsigned getNumberOfDerivatives(){ return 0; }
- void clearGrid();
- bool prepareForTasks();
+ void clearAverage();
+ void prepareForAveraging();
void compute( const unsigned& , MultiValue& ) const ;
};
@@ -230,7 +230,7 @@ ActionWithGrid(ao)
if( mycolv->isDensity() ) createGrid( "histogram", vstring );
else createGrid( "average", vstring );
// And finish the grid setup
- finishGridSetup();
+ setAveragingAction( mygrid, true );
// Enusre units for cube files are set correctly
if( !fractional ){
@@ -247,8 +247,7 @@ unsigned MultiColvarDensity::getNumberOfQuantities() const {
return directions.size() + 2;
}
-void MultiColvarDensity::clearGrid(){
- plumed_assert( mygrid->wasreset() );
+void MultiColvarDensity::clearAverage(){
std::vector<double> min(directions.size()), max(directions.size());
std::vector<std::string> gmin(directions.size()), gmax(directions.size());;
for(unsigned i=0;i<directions.size();++i){ min[i]=-0.5; max[i]=0.5; }
@@ -267,10 +266,11 @@ void MultiColvarDensity::clearGrid(){
}
}
for(unsigned i=0;i<directions.size();++i){ Tools::convert(min[i],gmin[i]); Tools::convert(max[i],gmax[i]); }
- mygrid->clear(); mygrid->setBounds( gmin, gmax, nbins, gspacing ); resizeFunctions();
+ ActionWithAveraging::clearAverage();
+ mygrid->setBounds( gmin, gmax, nbins, gspacing ); resizeFunctions();
}
-bool MultiColvarDensity::prepareForTasks(){
+void MultiColvarDensity::prepareForAveraging(){
for(unsigned i=0;i<directions.size();++i){
if( confined[i] ) continue;
std::string max; Tools::convert( 0.5*mycolv->getBox()(directions[i],directions[i]), max );
@@ -281,7 +281,7 @@ bool MultiColvarDensity::prepareForTasks(){
for(unsigned i=0;i<stash->getNumberOfStoredValues();++i) taskFlags[i]=1;
lockContributors();
// Retrieve the origin
- origin = getPosition(0); return true;
+ origin = getPosition(0);
}
void MultiColvarDensity::compute( const unsigned& current, MultiValue& myvals ) const {
diff --git a/src/vesselbase/ActionWithAveraging.cpp b/src/vesselbase/ActionWithAveraging.cpp
new file mode 100644
index 000000000..7a0b60340
--- /dev/null
+++ b/src/vesselbase/ActionWithAveraging.cpp
@@ -0,0 +1,126 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "ActionWithAveraging.h"
+#include "core/PlumedMain.h"
+#include "core/ActionSet.h"
+
+namespace PLMD {
+namespace vesselbase {
+
+void ActionWithAveraging::registerKeywords( Keywords& keys ){
+ Action::registerKeywords( keys ); ActionPilot::registerKeywords( keys );
+ ActionAtomistic::registerKeywords( keys ); ActionWithArguments::registerKeywords( keys ); ActionWithVessel::registerKeywords( keys );
+ keys.add("compulsory","STRIDE","1","the frequency with which the data should be collected and added to the quantity being averaged");
+ keys.add("compulsory","CLEAR","0","the frequency with which to clear all the accumulated data. The default value "
+ "of 0 implies that all the data will be used and that the grid will never be cleared");
+ keys.add("optional","LOGWEIGHTS","list of actions that calculates log weights that should be used to weight configurations when calculating averages");
+ keys.addFlag("UNORMALIZED",false,"output the unaveraged quantity/quantities.");
+}
+
+ActionWithAveraging::ActionWithAveraging( const ActionOptions& ao ):
+Action(ao),
+ActionPilot(ao),
+ActionAtomistic(ao),
+ActionWithArguments(ao),
+ActionWithVessel(ao),
+myaverage(NULL),
+useRunAllTasks(false),
+clearstride(0)
+{
+ if( keywords.exists("CLEAR") ){
+ parse("CLEAR",clearstride);
+ if( clearstride>0 ){
+ if( clearstride%getStride()!=0 ) error("CLEAR parameter must be a multiple of STRIDE");
+ log.printf(" clearing grid every %d steps \n",clearstride);
+ }
+ }
+ if( keywords.exists("LOGWEIGHTS") ){
+ std::vector<std::string> wwstr; parseVector("LOGWEIGHTS",wwstr);
+ log.printf(" reweighting using weights from ");
+ std::vector<Value*> arg( getArguments() );
+ for(unsigned i=0;i<wwstr.size();++i){
+ ActionWithValue* val = plumed.getActionSet().selectWithLabel<ActionWithValue*>(wwstr[i]);
+ if( !val ) error("could not find value named");
+ weights.push_back( val->copyOutput(val->getLabel()) );
+ arg.push_back( val->copyOutput(val->getLabel()) );
+ log.printf("%s ",wwstr[i].c_str() );
+ }
+ log.printf("\n"); requestArguments( arg );
+ }
+ if( keywords.exists("UNORMALIZED") ) parseFlag("UNORMALIZED",unormalised);
+}
+
+void ActionWithAveraging::setAveragingAction( AveragingVessel* av_vessel, const bool& usetasks ){
+ myaverage=av_vessel; addVessel( myaverage );
+ useRunAllTasks=usetasks; resizeFunctions();
+}
+
+void ActionWithAveraging::lockRequests(){
+ ActionAtomistic::lockRequests();
+ ActionWithArguments::lockRequests();
+}
+
+void ActionWithAveraging::unlockRequests(){
+ ActionAtomistic::unlockRequests();
+ ActionWithArguments::unlockRequests();
+}
+
+void ActionWithAveraging::calculateNumericalDerivatives(PLMD::ActionWithValue*){
+ error("not possible to compute numerical derivatives for this action");
+}
+
+void ActionWithAveraging::update(){
+ if( getStep()==0 || !onStep() ) return;
+ // Clear if it is time to reset
+ if( myaverage ){
+ if( myaverage->wasreset() ) clearAverage();
+ }
+ // Calculate the weight for all reweighting
+ if ( weights.size()>0 ){
+ double sum=0; for(unsigned i=0;i<weights.size();++i) sum+=weights[i]->get();
+ cweight = exp( sum );
+ } else {
+ cweight=1.0;
+ }
+ // Prepare to do the averaging
+ prepareForAveraging();
+ // Run all the tasks (if required
+ if( useRunAllTasks ) runAllTasks();
+ // This the averaging if it is not done using task list
+ else performOperations( true );
+ // Update the norm
+ if( myaverage ) myaverage->setNorm( cweight + myaverage->getNorm() );
+ // Finish the averaging
+ finishAveraging();
+ // By resetting here we are ensuring that the grid will be cleared at the start of the next step
+ if( myaverage ){
+ if( getStride()==0 || (clearstride>0 && getStep()%clearstride==0) ) myaverage->reset();
+ }
+}
+
+void ActionWithAveraging::clearAverage(){ plumed_assert( myaverage->wasreset() ); myaverage->clear(); }
+
+void ActionWithAveraging::performOperations( const bool& from_update ){ plumed_error(); }
+
+
+}
+}
diff --git a/src/vesselbase/ActionWithAveraging.h b/src/vesselbase/ActionWithAveraging.h
new file mode 100644
index 000000000..7642bcd17
--- /dev/null
+++ b/src/vesselbase/ActionWithAveraging.h
@@ -0,0 +1,94 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_vesselbase_ActionWithAveraging_h
+#define __PLUMED_vesselbase_ActionWithAveraging_h
+
+#include "core/ActionPilot.h"
+#include "core/ActionWithValue.h"
+#include "core/ActionAtomistic.h"
+#include "core/ActionWithArguments.h"
+#include "ActionWithVessel.h"
+#include "AveragingVessel.h"
+
+namespace PLMD {
+namespace vesselbase {
+
+/**
+\ingroup INHERIT
+This abstract base class should be used if you are writing some method that calculates an "average" from a set of
+trajectory frames. Notice that we use the word average very broadly here and state that even dimensionality
+reduction algorithms calculate an "average." In other words, what we mean by average is that the method is going
+to take in data from each trajectory frame and only calculate the final quantity once a certain amount of data has
+been collected.
+*/
+
+class ActionWithAveraging :
+ public ActionPilot,
+ public ActionAtomistic,
+ public ActionWithArguments,
+ public ActionWithVessel
+{
+friend class AveragingVessel;
+private:
+/// The vessel which is used to compute averages
+ AveragingVessel* myaverage;
+/// This ensures runAllTasks is used
+ bool useRunAllTasks;
+/// The frequency with which to clear the grid
+ unsigned clearstride;
+/// The weights we are going to use for reweighting
+ std::vector<Value*> weights;
+/// Are we accumulated the unormalized quantity
+ bool unormalised;
+protected:
+/// The current weight
+ double cweight;
+/// Set the averaging action
+ void setAveragingAction( AveragingVessel* av_vessel, const bool& usetasks );
+public:
+ static void registerKeywords( Keywords& keys );
+ explicit ActionWithAveraging( const ActionOptions& );
+ void lockRequests();
+ void unlockRequests();
+ void calculateNumericalDerivatives(PLMD::ActionWithValue*);
+ unsigned getNumberOfArguments() const ;
+ void calculate(){}
+ void apply(){}
+ void update();
+/// This does the clearing of the action
+ virtual void clearAverage();
+/// This is done before the averaging comences
+ virtual void prepareForAveraging(){}
+/// This does the averaging operation
+ virtual void performOperations( const bool& from_update );
+/// This is done once the averaging is finished
+ virtual void finishAveraging(){}
+};
+
+inline
+unsigned ActionWithAveraging::getNumberOfArguments() const {
+ return ActionWithArguments::getNumberOfArguments() - weights.size();
+}
+
+}
+}
+#endif
diff --git a/src/vesselbase/ActionWithVessel.cpp b/src/vesselbase/ActionWithVessel.cpp
index 876f89757..9bf4e72b7 100644
--- a/src/vesselbase/ActionWithVessel.cpp
+++ b/src/vesselbase/ActionWithVessel.cpp
@@ -294,8 +294,6 @@ void ActionWithVessel::runAllTasks(){
for(unsigned i=rank;i<nactive_tasks;i+=stride){
// Calculate the stuff in the loop for this action
performTask( indexOfTaskInFullList[i], partialTaskList[i], myvals );
- // Weight should be between zero and one
- plumed_dbg_assert( myvals.get(0)>=0 && myvals.get(0)<=1.0 );
// Check for conditions that allow us to just to skip the calculation
// the condition is that the weight of the contribution is low
diff --git a/src/vesselbase/AveragingVessel.cpp b/src/vesselbase/AveragingVessel.cpp
new file mode 100644
index 000000000..7d8933bd6
--- /dev/null
+++ b/src/vesselbase/AveragingVessel.cpp
@@ -0,0 +1,62 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "AveragingVessel.h"
+#include "ActionWithAveraging.h"
+
+namespace PLMD {
+namespace vesselbase {
+
+void AveragingVessel::registerKeywords( Keywords& keys ){
+ Vessel::registerKeywords( keys );
+}
+
+AveragingVessel::AveragingVessel( const vesselbase::VesselOptions& vo ):
+Vessel(vo),
+wascleared(true)
+{
+ ActionWithAveraging* myav = dynamic_cast<ActionWithAveraging*>( getAction() );
+ plumed_assert( myav ); unormalised = myav->unormalised;
+}
+
+void AveragingVessel::finish( const std::vector<double>& buffer ){
+ wascleared=false; for(unsigned i=1;i<data.size();++i) data[i]+=buffer[bufstart + i - 1];
+}
+
+bool AveragingVessel::wasreset() const {
+ return wascleared;
+}
+
+void AveragingVessel::clear(){
+ plumed_assert( wascleared ); data.assign( data.size(), 0.0 );
+}
+
+void AveragingVessel::reset(){
+ wascleared=true;
+}
+
+void AveragingVessel::setDataSize( const unsigned& size ){
+ if( data.size()!=(1+size) ) data.resize( 1+size, 0 );
+}
+
+}
+}
+
diff --git a/src/vesselbase/AveragingVessel.h b/src/vesselbase/AveragingVessel.h
new file mode 100644
index 000000000..c71427f93
--- /dev/null
+++ b/src/vesselbase/AveragingVessel.h
@@ -0,0 +1,101 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2012-2016 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#ifndef __PLUMED_vesselbase_AveragingVessel_h
+#define __PLUMED_vesselbase_AveragingVessel_h
+
+#include "Vessel.h"
+
+namespace PLMD {
+namespace vesselbase {
+
+class AveragingVessel : public Vessel {
+private:
+/// The grid was recently cleared and bounds can be set
+ bool wascleared;
+/// Are we outputting unormalised data
+ bool unormalised;
+/// The data that is being averaged
+ std::vector<double> data;
+protected:
+/// Set the size of the data vector
+ void setDataSize( const unsigned& size );
+/// Set an element of the data array
+ void setDataElement( const unsigned& myelem, const double& value );
+/// Add some value to an element of the data array
+ void addDataElement( const unsigned& myelem, const double& value );
+/// Get the value of one of the data element
+ double getDataElement( const unsigned& myelem ) const ;
+/// Are we averaging the data
+ bool noAverage() const { return unormalised; }
+public:
+/// keywords
+ static void registerKeywords( Keywords& keys );
+/// Constructor
+ explicit AveragingVessel( const vesselbase::VesselOptions& );
+/// Copy data from an accumulated buffer into the grid
+ virtual void finish( const std::vector<double>& );
+/// Was the grid cleared on the last step
+ bool wasreset() const ;
+/// Clear all the data stored on the grid
+ virtual void clear();
+/// Reset the grid so that it is cleared at start of next time it is calculated
+ virtual void reset();
+/// Functions for dealing with normalisation constant
+ void setNorm( const double& snorm );
+ double getNorm() const ;
+ bool applyForce( std::vector<double>& forces ){ return false; }
+};
+
+inline
+void AveragingVessel::setDataElement( const unsigned& myelem, const double& value ){
+ plumed_dbg_assert( myelem<1+data.size() );
+ wascleared=false; data[1+myelem]=value;
+}
+
+inline
+void AveragingVessel::addDataElement( const unsigned& myelem, const double& value ){
+ plumed_dbg_assert( myelem<1+data.size() );
+ wascleared=false; data[1+myelem]+=value;
+}
+
+inline
+double AveragingVessel::getDataElement( const unsigned& myelem ) const {
+ plumed_dbg_assert( myelem<data.size()-1 );
+ if( unormalised ) return data[1+myelem];
+ return data[1+myelem] / data[0];
+}
+
+inline
+void AveragingVessel::setNorm( const double& snorm ){
+ plumed_dbg_assert( data.size()>0 );
+ data[0]=snorm;
+}
+
+inline
+double AveragingVessel::getNorm() const {
+ plumed_dbg_assert( data.size()>0 );
+ return data[0];
+}
+
+}
+}
+#endif
--
GitLab