From 190093f94b3dd49a56646acc8cfe5bed601c7737 Mon Sep 17 00:00:00 2001 From: Carlo Camilloni Date: Thu, 19 Nov 2020 07:03:01 +0100 Subject: [PATCH] WIP: MetaD NEIGHBOR list option to sum gaussians (#639) * MetaD: this commits add the neighbor list option to sum hills this can be usefull with many cvs where grids are either too memory consuming and/or too slow. I was playing with this idea with OPES and it actually works pretty well. The neighbor list update at the moment is only openMP parallelised. this commits add also the openMP parallelisation to the cycle over the gaussians * MetaD: some fixes for the OMP implementation * MetaD NL, not use mpi for the neighbor list update * MetaD: split getBiasDer and not Der * MetaD: more cleaning and NEIGHBOR restart * MetaD: added missing components in doc * MetaD: some reorganisation and fix * MetaD: added control for too few hills and openMP * MetaD: take also into account mpi threads to avoid the use of too many threads for few hills * increased the number of gaussians in two regtest to cover better the parallelism * Fixed grids * Metad: cleaning * MetaD: performanaces and cleaning * more std * added back std::time * Metad: fix thread safety without loss of performances * MetaD: some cleaning in the initialisation * METAD: regtest for CALC_WORK and NLIST * Metad: added checks that RESTART and WALKERS_MPI behave consistently among the walkers * Metad: cleaning * Metad: updated documentation and added NLIST_PARAMETERS to modify the neighbor list rules * Metad: slightly improved performances of GRID and more docs * MetaD: moved a line in nlistUpdate * Metad: spelling * MetaD: in case of WALKERS_MPI the root replica will set the folder from which reading GRID or HILLS files upon RESTART * doc * Metad - restart stops if hills not found - broadcasting avoided if absolute path provided * Metad: restart with hills give error if not already restarted from grid * fixes for new METAD restart checks * fixed regteast * Metad: removed unused headers added two missing std --- CHANGES/v2.8.md | 8 + regtest/basic/rt-mpi7/gridx.0.reference | 1 - regtest/basic/rt-mpi7/plumed.dat | 2 + regtest/basic/rt-mpi7b/grid.0.reference | 2 +- regtest/basic/rt-mpi7b/grid.1.reference | 2 +- regtest/basic/rt-mpi7b/grid.2.reference | 2 +- regtest/basic/rt-mpi7b/gridx.0.reference | 11 +- regtest/basic/rt-mpi7b/gridx.1.reference | 10 +- regtest/basic/rt-mpi7b/gridx.2.reference | 10 +- regtest/basic/rt-mpi7b/plumed.dat | 2 - regtest/basic/rt-mpi9/plumed.dat | 2 + regtest/basic/rt10-mpi/plumed.dat | 2 +- regtest/basic/rt10-restart/COLVAR.reference | 12 +- regtest/basic/rt10-restart/H2 | 13 + regtest/basic/rt10-restart/plumed.dat | 6 +- regtest/basic/rt10/COLVAR.reference | 6 +- regtest/basic/rt10/HILLS.reference | 2 - regtest/basic/rt10/forces.reference | 18 +- regtest/basic/rt10/plumed.dat | 3 +- regtest/basic/rt81/plumed.dat | 2 +- .../rt-funnel-mwalkers/COLVAR.0.reference | 62 +- .../rt-funnel-mwalkers/COLVAR.1.reference | 62 +- regtest/funnel/rt-funnel/COLVAR.reference | 62 +- regtest/funnel/rt-sphere/COLVAR.reference | 62 +- src/bias/MetaD.cpp | 1142 ++++++++++------- 25 files changed, 897 insertions(+), 609 deletions(-) create mode 100644 regtest/basic/rt10-restart/H2 diff --git a/CHANGES/v2.8.md b/CHANGES/v2.8.md index 8a20259971..73272fe916 100644 --- a/CHANGES/v2.8.md +++ b/CHANGES/v2.8.md @@ -4,6 +4,14 @@ This page contains changes that will end up in 2.8 +- Changes leading to differences with previous versions + - in \ref METAD if possible the root walker in WALKERS_MPI will set the folder from which reading the GRID/HILLS file upon restart + - in \ref METAD work is not calculated by default anymore, if needed it can be obtained using CALC_WORK + +- Other improvements + - in \ref METAD a new keyword NLIST has been added to use a neighbor list for bias evaluation, this should be faster than grids with many CVs + - in \ref METAD there are more checks that a restart of WALKERS_MPI is working consistently among walkers + - Changes in the OPES module - new action \ref OPES_EXPANDED - various new actions of type \ref EXPANSION_CV to be used with \ref OPES_EXPANDED diff --git a/regtest/basic/rt-mpi7/gridx.0.reference b/regtest/basic/rt-mpi7/gridx.0.reference index 624de6d661..c317f66b74 100644 --- a/regtest/basic/rt-mpi7/gridx.0.reference +++ b/regtest/basic/rt-mpi7/gridx.0.reference @@ -1,4 +1,3 @@ -# this is to test restart #! FIELDS d1 d2 meta2.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 diff --git a/regtest/basic/rt-mpi7/plumed.dat b/regtest/basic/rt-mpi7/plumed.dat index 20044e1cab..7a5961554a 100644 --- a/regtest/basic/rt-mpi7/plumed.dat +++ b/regtest/basic/rt-mpi7/plumed.dat @@ -9,6 +9,7 @@ METAD ... GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4 GRID_WFILE=grid GRID_WSTRIDE=2 LABEL=meta1 + RESTART=NO ... METAD ... @@ -18,6 +19,7 @@ METAD ... GRID_WFILE=gridx GRID_WSTRIDE=2 STORE_GRIDS FILE=HILLSX LABEL=meta2 + RESTART=NO ... METAD ... diff --git a/regtest/basic/rt-mpi7b/grid.0.reference b/regtest/basic/rt-mpi7b/grid.0.reference index eba63d858b..85ce651b49 100644 --- a/regtest/basic/rt-mpi7b/grid.0.reference +++ b/regtest/basic/rt-mpi7b/grid.0.reference @@ -1,4 +1,4 @@ -#! FIELDS d1 d2 @3.bias der_d1 der_d2 +#! FIELDS d1 d2 @2.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/grid.1.reference b/regtest/basic/rt-mpi7b/grid.1.reference index f68da103b2..b166215c90 100644 --- a/regtest/basic/rt-mpi7b/grid.1.reference +++ b/regtest/basic/rt-mpi7b/grid.1.reference @@ -1,4 +1,4 @@ -#! FIELDS d1 d2 @3.bias der_d1 der_d2 +#! FIELDS d1 d2 @2.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/grid.2.reference b/regtest/basic/rt-mpi7b/grid.2.reference index eb0d399a96..c0754dc6ff 100644 --- a/regtest/basic/rt-mpi7b/grid.2.reference +++ b/regtest/basic/rt-mpi7b/grid.2.reference @@ -1,4 +1,4 @@ -#! FIELDS d1 d2 @3.bias der_d1 der_d2 +#! FIELDS d1 d2 @2.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/gridx.0.reference b/regtest/basic/rt-mpi7b/gridx.0.reference index 1dd64c8f88..79880ff4a1 100644 --- a/regtest/basic/rt-mpi7b/gridx.0.reference +++ b/regtest/basic/rt-mpi7b/gridx.0.reference @@ -1,5 +1,4 @@ -# this is to test restart -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -37,7 +36,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -75,7 +74,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -113,7 +112,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -151,7 +150,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/gridx.1.reference b/regtest/basic/rt-mpi7b/gridx.1.reference index f6f7385b4c..ff845b1bfe 100644 --- a/regtest/basic/rt-mpi7b/gridx.1.reference +++ b/regtest/basic/rt-mpi7b/gridx.1.reference @@ -1,4 +1,4 @@ -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -36,7 +36,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -74,7 +74,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -112,7 +112,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -150,7 +150,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/gridx.2.reference b/regtest/basic/rt-mpi7b/gridx.2.reference index aaccbbd58c..1cdc3e454e 100644 --- a/regtest/basic/rt-mpi7b/gridx.2.reference +++ b/regtest/basic/rt-mpi7b/gridx.2.reference @@ -1,4 +1,4 @@ -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -36,7 +36,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -74,7 +74,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -112,7 +112,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 @@ -150,7 +150,7 @@ 1.000000000 2.000000000 0.000000000 0.000000000 0.000000000 1.500000000 2.000000000 0.000000000 0.000000000 0.000000000 2.000000000 2.000000000 0.000000000 0.000000000 0.000000000 -#! FIELDS d1 d2 @4.bias der_d1 der_d2 +#! FIELDS d1 d2 @3.bias der_d1 der_d2 #! SET min_d1 0 #! SET max_d1 2 #! SET nbins_d1 5 diff --git a/regtest/basic/rt-mpi7b/plumed.dat b/regtest/basic/rt-mpi7b/plumed.dat index 36bc295a9d..710f15055a 100644 --- a/regtest/basic/rt-mpi7b/plumed.dat +++ b/regtest/basic/rt-mpi7b/plumed.dat @@ -1,5 +1,3 @@ -RESTART - d1: DISTANCE ATOMS=1,2 d2: DISTANCE ATOMS=3,4 diff --git a/regtest/basic/rt-mpi9/plumed.dat b/regtest/basic/rt-mpi9/plumed.dat index af5bbaacca..5a1e2dd1f5 100644 --- a/regtest/basic/rt-mpi9/plumed.dat +++ b/regtest/basic/rt-mpi9/plumed.dat @@ -8,6 +8,7 @@ METAD ... FMT=%6.2f GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4 GRID_WFILE=grid GRID_WSTRIDE=2 + RESTART=NO ... METAD ... @@ -16,6 +17,7 @@ METAD ... GRID_MIN=0,0 GRID_MAX=2,2 GRID_BIN=4,4 GRID_WFILE=gridx GRID_WSTRIDE=2 STORE_GRIDS FILE=HILLSX + RESTART=NO ... diff --git a/regtest/basic/rt10-mpi/plumed.dat b/regtest/basic/rt10-mpi/plumed.dat index ebde2b3ad5..be9254f3e6 100644 --- a/regtest/basic/rt10-mpi/plumed.dat +++ b/regtest/basic/rt10-mpi/plumed.dat @@ -1,7 +1,7 @@ mu1: DISTANCE ATOMS=1,10 vol: VOLUME -METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=10 LABEL=md FMT=%14.6f +METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=2 LABEL=md FMT=%14.6f PRINT ... STRIDE=1 diff --git a/regtest/basic/rt10-restart/COLVAR.reference b/regtest/basic/rt10-restart/COLVAR.reference index f430a3170f..b06b9ef451 100644 --- a/regtest/basic/rt10-restart/COLVAR.reference +++ b/regtest/basic/rt10-restart/COLVAR.reference @@ -4,9 +4,9 @@ 0.100000 1.098 127.933 1.759 0.150000 1.080 127.933 2.578 0.200000 1.087 127.933 3.651 -#! FIELDS time mu1 vol md.bias - 0.000000 1.163 127.933 4.223 - 0.050000 1.131 127.933 4.688 - 0.100000 1.098 127.933 5.686 - 0.150000 1.080 127.933 6.441 - 0.200000 1.087 127.933 7.552 +#! FIELDS time mu1 vol md.work md.bias md1.bias md2.bias md2.nlker md2.nlsteps + 0.000000 1.163 127.933 0.000 4.223 0.000 0.000 9.000 1.000 + 0.050000 1.131 127.933 0.000 4.688 0.000 0.000 9.000 1.000 + 0.100000 1.098 127.933 1.000 5.686 0.000 0.000 10.000 2.000 + 0.150000 1.080 127.933 2.000 6.441 0.000 0.000 11.000 1.000 + 0.200000 1.087 127.933 3.000 7.552 0.000 0.000 12.000 1.000 diff --git a/regtest/basic/rt10-restart/H2 b/regtest/basic/rt10-restart/H2 new file mode 100644 index 0000000000..cc37a62c8a --- /dev/null +++ b/regtest/basic/rt10-restart/H2 @@ -0,0 +1,13 @@ +#! FIELDS time mu1 vol sigma_mu1 sigma_vol height biasf +#! SET multivariate false + 0 1.162646040831079 127.932640011072 0.1 0.2 0 1 + 0.05 1.130546273059004 127.932640011072 0.1 0.2 0 1 + 0.1 1.097928292824707 127.932640011072 0.1 0.2 0 1 + 0.15 1.080244153391634 127.932640011072 0.1 0.2 0 1 + 0.2 1.086854650075657 127.932640011072 0.1 0.2 0 1 +#! FIELDS time mu1 vol sigma_mu1 sigma_vol height biasf +#! SET multivariate false + 0.050000 1.130546 127.932640 0.100000 0.200000 0.000000 1.000000 + 0.100000 1.097928 127.932640 0.100000 0.200000 0.000000 1.000000 + 0.150000 1.080244 127.932640 0.100000 0.200000 0.000000 1.000000 + 0.200000 1.086855 127.932640 0.100000 0.200000 0.000000 1.000000 diff --git a/regtest/basic/rt10-restart/plumed.dat b/regtest/basic/rt10-restart/plumed.dat index e8d2428b7d..c6e0115bbd 100644 --- a/regtest/basic/rt10-restart/plumed.dat +++ b/regtest/basic/rt10-restart/plumed.dat @@ -2,15 +2,17 @@ RESTART mu1: DISTANCE ATOMS=1,10 vol: VOLUME -METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=10 LABEL=md FMT=%14.6f +METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=10 LABEL=md FMT=%14.6f CALC_WORK # this is to check RESTART. # since here restart is disabled, the hills in H1 should not be read METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=0.0 PACE=10 LABEL=md1 FMT=%14.6f FILE=H1 RESTART=NO +METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=0.0 PACE=10 LABEL=md2 FMT=%14.6f NLIST FILE=H2 + PRINT ... STRIDE=1 - ARG=mu1,vol,md.bias + ARG=mu1,vol,md.work,md.bias,md1.bias,md2.* FILE=COLVAR FMT=%6.3f ... PRINT diff --git a/regtest/basic/rt10/COLVAR.reference b/regtest/basic/rt10/COLVAR.reference index 7504a6077a..347b7eb886 100644 --- a/regtest/basic/rt10/COLVAR.reference +++ b/regtest/basic/rt10/COLVAR.reference @@ -1,6 +1,6 @@ #! FIELDS time mu1 vol md.bias 0.000000 1.163 127.933 0.000 0.050000 1.131 127.933 0.000 - 0.100000 1.098 127.933 0.948 - 0.150000 1.080 127.933 1.866 - 0.200000 1.087 127.933 2.901 + 0.100000 1.098 127.933 0.000 + 0.150000 1.080 127.933 0.984 + 0.200000 1.087 127.933 0.994 diff --git a/regtest/basic/rt10/HILLS.reference b/regtest/basic/rt10/HILLS.reference index 7564525a0c..b6ddf9f9d5 100644 --- a/regtest/basic/rt10/HILLS.reference +++ b/regtest/basic/rt10/HILLS.reference @@ -1,7 +1,5 @@ #! FIELDS time mu1 vol sigma_mu1 sigma_vol height biasf #! SET multivariate false #! SET kerneltype gaussian - 0.050000 1.130546 127.932640 0.100000 0.200000 1.000000 -1.000000 0.100000 1.097928 127.932640 0.100000 0.200000 1.000000 -1.000000 - 0.150000 1.080244 127.932640 0.100000 0.200000 1.000000 -1.000000 0.200000 1.086855 127.932640 0.100000 0.200000 1.000000 -1.000000 diff --git a/regtest/basic/rt10/forces.reference b/regtest/basic/rt10/forces.reference index 0f6df3cbbd..0b8cb3bdfb 100644 --- a/regtest/basic/rt10/forces.reference +++ b/regtest/basic/rt10/forces.reference @@ -219,8 +219,9 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 108 - 1.498973 0.002364 1.894352 -X 2.054883 -0.081609 -2.310045 + 0.000000 0.000000 0.000000 +X 0.000000 0.000000 0.000000 +X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 @@ -229,7 +230,6 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 -X -2.054883 0.081609 2.310045 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 @@ -329,8 +329,8 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 108 - 2.617465 0.034167 4.017150 -X 3.867602 -0.441880 -4.791378 + 0.738158 0.009635 1.132887 +X 1.090712 -0.124616 -1.351229 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 @@ -339,7 +339,7 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 -X -3.867602 0.441880 4.791378 +X -1.090712 0.124616 1.351229 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 @@ -439,8 +439,8 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 108 - 1.570522 0.052663 3.172461 -X 2.525074 -0.462386 -3.588808 + 0.391738 0.013136 0.791313 +X 0.629834 -0.115334 -0.895163 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 @@ -449,7 +449,7 @@ X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 -X -2.525074 0.462386 3.588808 +X -0.629834 0.115334 0.895163 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 X 0.000000 0.000000 0.000000 diff --git a/regtest/basic/rt10/plumed.dat b/regtest/basic/rt10/plumed.dat index 3edaf19466..acf26746f7 100644 --- a/regtest/basic/rt10/plumed.dat +++ b/regtest/basic/rt10/plumed.dat @@ -1,8 +1,7 @@ mu1: DISTANCE ATOMS=1,10 vol: VOLUME -METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=5*2 LABEL=md FMT=%14.6f - +METAD ARG=mu1,vol SIGMA=0.1,0.2 HEIGHT=1.0 PACE=2*2 LABEL=md FMT=%14.6f PRINT ... STRIDE=1 ARG=mu1,vol,md.bias diff --git a/regtest/basic/rt81/plumed.dat b/regtest/basic/rt81/plumed.dat index 69faa3ade3..2861aecc83 100644 --- a/regtest/basic/rt81/plumed.dat +++ b/regtest/basic/rt81/plumed.dat @@ -1,7 +1,7 @@ d: DISTANCE ATOMS=1,10 c: COORDINATION GROUPA=1-108 GROUPB=1-108 R_0=0.5 -fg: METAD ARG=c,d SIGMA=0.1,0.2 HEIGHT=5.0 PACE=1 FMT=%14.6f WALKERS_MPI FILE=HILLS_multi FLYING_GAUSSIAN RESTART=YES +fg: METAD ARG=c,d SIGMA=0.1,0.2 HEIGHT=5.0 PACE=1 FMT=%14.6f WALKERS_MPI FILE=HILLS_multi FLYING_GAUSSIAN PRINT ARG=fg.bias FILE=colvar diff --git a/regtest/funnel/rt-funnel-mwalkers/COLVAR.0.reference b/regtest/funnel/rt-funnel-mwalkers/COLVAR.0.reference index d11c2d0d30..38e39cc785 100644 --- a/regtest/funnel/rt-funnel-mwalkers/COLVAR.0.reference +++ b/regtest/funnel/rt-funnel-mwalkers/COLVAR.0.reference @@ -1,31 +1,31 @@ -#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias metad.work - 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 0.000000 - 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 0.000000 - 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 0.000000 - 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 0.000000 - 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 0.000000 - 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 0.000000 - 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 0.000000 - 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 0.000000 - 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 0.000000 - 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 0.000000 - 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 0.000000 - 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 0.000000 - 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 0.000000 - 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 0.000000 - 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 0.000000 - 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 0.000000 - 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 0.000000 - 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 0.000000 - 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 0.000000 - 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 0.000000 - 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 0.000000 - 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 0.000000 - 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 0.000000 - 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 0.000000 - 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 0.000000 - 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 0.000000 - 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 0.000000 - 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 0.000000 - 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 0.000000 - 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 0.000000 +#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias + 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 + 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 + 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 + 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 + 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 + 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 + 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 + 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 + 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 + 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 + 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 + 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 + 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 + 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 + 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 + 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 + 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 + 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 + 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 + 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 + 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 + 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 + 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 + 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 + 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 + 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 + 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 + 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 + 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 + 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 diff --git a/regtest/funnel/rt-funnel-mwalkers/COLVAR.1.reference b/regtest/funnel/rt-funnel-mwalkers/COLVAR.1.reference index d11c2d0d30..38e39cc785 100644 --- a/regtest/funnel/rt-funnel-mwalkers/COLVAR.1.reference +++ b/regtest/funnel/rt-funnel-mwalkers/COLVAR.1.reference @@ -1,31 +1,31 @@ -#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias metad.work - 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 0.000000 - 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 0.000000 - 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 0.000000 - 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 0.000000 - 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 0.000000 - 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 0.000000 - 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 0.000000 - 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 0.000000 - 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 0.000000 - 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 0.000000 - 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 0.000000 - 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 0.000000 - 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 0.000000 - 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 0.000000 - 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 0.000000 - 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 0.000000 - 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 0.000000 - 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 0.000000 - 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 0.000000 - 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 0.000000 - 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 0.000000 - 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 0.000000 - 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 0.000000 - 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 0.000000 - 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 0.000000 - 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 0.000000 - 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 0.000000 - 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 0.000000 - 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 0.000000 - 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 0.000000 +#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias + 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 + 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 + 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 + 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 + 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 + 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 + 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 + 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 + 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 + 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 + 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 + 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 + 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 + 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 + 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 + 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 + 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 + 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 + 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 + 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 + 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 + 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 + 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 + 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 + 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 + 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 + 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 + 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 + 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 + 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 diff --git a/regtest/funnel/rt-funnel/COLVAR.reference b/regtest/funnel/rt-funnel/COLVAR.reference index d11c2d0d30..38e39cc785 100644 --- a/regtest/funnel/rt-funnel/COLVAR.reference +++ b/regtest/funnel/rt-funnel/COLVAR.reference @@ -1,31 +1,31 @@ -#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias metad.work - 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 0.000000 - 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 0.000000 - 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 0.000000 - 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 0.000000 - 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 0.000000 - 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 0.000000 - 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 0.000000 - 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 0.000000 - 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 0.000000 - 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 0.000000 - 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 0.000000 - 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 0.000000 - 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 0.000000 - 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 0.000000 - 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 0.000000 - 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 0.000000 - 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 0.000000 - 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 0.000000 - 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 0.000000 - 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 0.000000 - 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 0.000000 - 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 0.000000 - 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 0.000000 - 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 0.000000 - 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 0.000000 - 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 0.000000 - 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 0.000000 - 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 0.000000 - 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 0.000000 - 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 0.000000 +#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias + 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 + 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 + 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 + 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 + 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 + 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 + 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 + 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 + 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 + 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 + 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 + 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 + 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 + 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 + 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 + 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 + 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 + 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 + 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 + 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 + 20.000000 1.925427 1.237743 0.657797 767.112807 0.000000 + 21.000000 2.033292 1.295877 0.656577 1041.494304 0.000000 + 22.000000 1.444796 0.849947 0.838569 420.826701 0.000000 + 23.000000 3.514052 2.780271 0.599456 4350.448397 0.000000 + 24.000000 1.975710 1.353646 0.685574 1658.456302 0.000000 + 25.000000 2.470418 1.763036 0.542308 3113.860056 0.000000 + 26.000000 3.169386 2.503920 0.676148 5793.739791 0.000000 + 27.000000 2.224747 1.762360 0.622640 4427.116046 0.000000 + 28.000000 2.660713 1.995432 0.681466 5894.334523 0.000000 + 29.000000 3.248468 2.525448 0.625886 4841.111588 0.000000 diff --git a/regtest/funnel/rt-sphere/COLVAR.reference b/regtest/funnel/rt-sphere/COLVAR.reference index a35f0a5363..7e5793f81a 100644 --- a/regtest/funnel/rt-sphere/COLVAR.reference +++ b/regtest/funnel/rt-sphere/COLVAR.reference @@ -1,31 +1,31 @@ -#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias metad.work - 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 0.000000 - 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 0.000000 - 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 0.000000 - 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 0.000000 - 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 0.000000 - 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 0.000000 - 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 0.000000 - 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 0.000000 - 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 0.000000 - 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 0.000000 - 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 0.000000 - 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 0.000000 - 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 0.000000 - 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 0.000000 - 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 0.000000 - 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 0.000000 - 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 0.000000 - 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 0.000000 - 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 0.000000 - 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 0.000000 - 20.000000 1.925427 1.237743 0.657797 0.000000 0.000000 0.000000 - 21.000000 2.033292 1.295877 0.656577 0.000000 0.000000 0.000000 - 22.000000 1.444796 0.849947 0.838569 0.000000 0.000000 0.000000 - 23.000000 3.514052 2.780271 0.599456 5814.922759 0.000000 0.000000 - 24.000000 1.975710 1.353646 0.685574 0.000000 0.000000 0.000000 - 25.000000 2.470418 1.763036 0.542308 38.168391 0.000000 0.000000 - 26.000000 3.169386 2.503920 0.676148 7850.367690 0.000000 0.000000 - 27.000000 2.224747 1.762360 0.622640 87.067878 0.000000 0.000000 - 28.000000 2.660713 1.995432 0.681466 1677.446418 0.000000 0.000000 - 29.000000 3.248468 2.525448 0.625886 6472.223271 0.000000 0.000000 +#! FIELDS time d1 fps.lp fps.ld funnel.bias metad.bias + 0.000000 0.396158 -0.009720 0.626527 0.000000 0.000000 + 1.000000 0.391878 -0.002473 0.559116 0.000000 0.000000 + 2.000000 0.805812 0.142132 0.533798 0.000000 0.000000 + 3.000000 0.849127 0.200841 0.548150 0.000000 0.000000 + 4.000000 0.730758 0.110941 0.572736 0.000000 0.000000 + 5.000000 1.138845 0.236602 0.584752 0.000000 0.000000 + 6.000000 1.258723 0.320768 0.525443 0.000000 0.000000 + 7.000000 1.162686 0.265523 0.581793 0.000000 0.000000 + 8.000000 0.903258 0.075093 0.572080 0.000000 0.000000 + 9.000000 0.858768 0.018716 0.626715 0.000000 0.000000 + 10.000000 0.849164 0.047239 0.604750 0.000000 0.000000 + 11.000000 1.107866 0.104703 0.613068 0.000000 0.000000 + 12.000000 1.097381 0.187387 0.533305 0.000000 0.000000 + 13.000000 1.291115 0.315166 0.549732 0.000000 0.000000 + 14.000000 1.298709 0.388939 0.529146 0.000000 0.000000 + 15.000000 1.302892 0.331084 0.528228 0.000000 0.000000 + 16.000000 1.354794 0.407920 0.617995 0.000000 0.000000 + 17.000000 1.292765 0.310389 0.454396 0.000000 0.000000 + 18.000000 1.431399 0.404609 0.301176 0.000000 0.000000 + 19.000000 1.394466 0.404893 0.244238 0.000000 0.000000 + 20.000000 1.925427 1.237743 0.657797 0.000000 0.000000 + 21.000000 2.033292 1.295877 0.656577 0.000000 0.000000 + 22.000000 1.444796 0.849947 0.838569 0.000000 0.000000 + 23.000000 3.514052 2.780271 0.599456 5814.922759 0.000000 + 24.000000 1.975710 1.353646 0.685574 0.000000 0.000000 + 25.000000 2.470418 1.763036 0.542308 38.168391 0.000000 + 26.000000 3.169386 2.503920 0.676148 7850.367690 0.000000 + 27.000000 2.224747 1.762360 0.622640 87.067878 0.000000 + 28.000000 2.660713 1.995432 0.681466 1677.446418 0.000000 + 29.000000 3.248468 2.525448 0.625886 6472.223271 0.000000 diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp index e8a08752cc..ee13cdd9e9 100644 --- a/src/bias/MetaD.cpp +++ b/src/bias/MetaD.cpp @@ -22,26 +22,23 @@ #include "Bias.h" #include "ActionRegister.h" #include "core/ActionSet.h" -#include "tools/Grid.h" #include "core/PlumedMain.h" #include "core/Atoms.h" -#include "tools/Exception.h" #include "core/FlexibleBin.h" +#include "tools/Exception.h" +#include "tools/Grid.h" #include "tools/Matrix.h" +#include "tools/OpenMP.h" #include "tools/Random.h" -#include -#include #include "tools/File.h" -#include -#include #include -#include +#include +#if defined(__PLUMED_HAS_GETCWD) +#include +#endif #define DP2CUTOFF 6.25 -using namespace std; - - namespace PLMD { namespace bias { @@ -78,18 +75,27 @@ calculation increases with the length of the simulation as one has to, at every the values of a larger and larger number of Gaussian kernels. To avoid this issue you can store the bias on a grid. This approach is similar to that proposed in \cite babi08jcp but has the advantage that the grid spacing is independent on the Gaussian width. -Notice that you should -provide either the number of bins for every collective variable (GRID_BIN) or -the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use -the most conservative choice (highest number of bins) for each dimension. +Notice that you should provide the grid boundaries (GRID_MIN and GRID_MAX) and either the number of bins +for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING). +In case you provide both PLUMED will use the most conservative choice (highest number of bins) for each dimension. In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING) -and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing. -This default choice should be reasonable for most applications. +PLUMED will use 1/5 of the Gaussian width (SIGMA) as grid spacing if the width is fixed or 1/5 of the minimum +Gaussian width (SIGMA_MIN) if the width is variable. This default choice should be reasonable for most applications. + +Alternatively to the use of grids, it is possible to use a neighbor list to decrease the cost of evaluating the bias, +this can be enabled using NLIST. NLIST can be beneficial with more than 2 collective variables, where GRID becomes +expensive and memory consuming. The neighbor list will be updated everytime the CVs go farther than a cut-off value +from the position they were at last neighbor list update. Gaussians are added to the neigbhor list if their center +is within 6.*DP2CUTOFF*sigma*sigma. While the list is updated if the CVs are farther from the center than 0.5 of the +standard deviation of the Gaussian center distribution of the list. These parameters (6 and 0.5) can be modified using +NLIST_PARAMETERS. Note that the use of neighbor list does not provide the exact bias. Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read it using GRID_RFILE. +The work performed by the METAD bias can be calculated using CALC_WORK, note that this is expensive when not using grids. + Another option that is available in plumed is well-tempered metadynamics \cite Barducci:2008. In this variant of metadynamics the heights of the Gaussian hills are scaled at each step so the bias is now given by: @@ -364,10 +370,6 @@ METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0 \endplumedfile The number of elements in the RECT array should be equal to the number of replicas. - - - - */ //+ENDPLUMEDOC @@ -375,15 +377,18 @@ class MetaD : public Bias { private: struct Gaussian { - vector center; - vector sigma; - double height; bool multivariate; // this is required to discriminate the one dimensional case - vector invsigma; - Gaussian(const vector & center,const vector & sigma,double height, bool multivariate ): - center(center),sigma(sigma),height(height),multivariate(multivariate),invsigma(sigma) { + double height; + std::vector center; + std::vector sigma; + std::vector invsigma; + Gaussian(const bool m, const double h, const std::vector& c, const std::vector& s): + multivariate(m),height(h),center(c),sigma(s),invsigma(s) { // to avoid troubles from zero element in flexible hills - for(unsigned i=0; i1.e-20) invsigma[i]=1.0/invsigma[i] ; else invsigma[i]=0.0; + for(unsigned i=0; i1.e-20) invsigma[i]=1.0/invsigma[i] ; + else invsigma[i]=0.0; + } } }; struct TemperingSpecs { @@ -397,82 +402,105 @@ class MetaD : public Bias { is_active(is_active), name_stem(name_stem), name(name), biasf(biasf), threshold(threshold), alpha(alpha) {} }; - vector sigma0_; - vector sigma0min_; - vector sigma0max_; - vector hills_; - OFile hillsOfile_; - OFile gridfile_; - std::unique_ptr BiasGrid_; - bool storeOldGrids_; - int wgridstride_; - bool grid_; - double height0_; - double biasf_; - static const size_t n_tempering_options_ = 1; - static const string tempering_names_[1][2]; - double dampfactor_; - struct TemperingSpecs tt_specs_; - std::string targetfilename_; - std::unique_ptr TargetGrid_; + // general setup double kbt_; int stride_; + bool calc_work_; + // well-tempered MetaD bool welltemp_; - // - int current_stride; - bool freq_adaptive_; - int fa_update_frequency_; - int fa_max_stride_; - double fa_min_acceleration_; - // - std::unique_ptr dp_; + double biasf_; + // output files format + std::string fmt_; + // first step? + bool isFirstStep_; + // Gaussian starting parameters + double height0_; + std::vector sigma0_; + std::vector sigma0min_; + std::vector sigma0max_; + // Gaussians + std::vector hills_; + std::unique_ptr flexbin_; int adaptive_; - std::unique_ptr flexbin; + OFile hillsOfile_; + std::vector> ifiles_; + std::vector ifilesnames_; + // Grids + bool grid_; + std::unique_ptr BiasGrid_; + OFile gridfile_; + bool storeOldGrids_; + int wgridstride_; + // multiple walkers int mw_n_; - string mw_dir_; + std::string mw_dir_; int mw_id_; int mw_rstride_; - bool walkers_mpi; + bool walkers_mpi_; unsigned mpi_nw_; unsigned mpi_mw_; - bool flying; - bool acceleration; - double acc; + // flying gaussians + bool flying_; + // kinetics from metadynamics + bool acceleration_; + double acc_; double acc_restart_mean_; + // transition-tempering metadynamics bool calc_max_bias_; double max_bias_; bool calc_transition_bias_; double transition_bias_; - vector > transitionwells_; - vector> ifiles; - vector ifilesnames; + std::vector > transitionwells_; + static const size_t n_tempering_options_ = 1; + static const std::string tempering_names_[1][2]; + double dampfactor_; + struct TemperingSpecs tt_specs_; + std::string targetfilename_; + std::unique_ptr TargetGrid_; + // frequency adaptive metadynamics + int current_stride_; + bool freq_adaptive_; + int fa_update_frequency_; + int fa_max_stride_; + double fa_min_acceleration_; + // intervals double uppI_; double lowI_; bool doInt_; - bool isFirstStep; + // reweighting bool calc_rct_; double reweight_factor_; unsigned rct_ustride_; + // work double work_; - long int last_step_warn_grid; - - static void registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys); + // neighbour list stuff + bool nlist_; + bool nlist_update_; + unsigned nlist_steps_; + std::array nlist_param_; + std::vector nlist_hills_; + std::vector nlist_center_; + std::vector nlist_dev2_; + + static void registerTemperingKeywords(const std::string &name_stem, const std::string &name, Keywords &keys); void readTemperingSpecs(TemperingSpecs &t_specs); void logTemperingSpecs(const TemperingSpecs &t_specs); void readGaussians(IFile*); void writeGaussian(const Gaussian&,OFile&); void addGaussian(const Gaussian&); - double getHeight(const vector&); + double getHeight(const std::vector&); void temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias); - double getBiasAndDerivatives(const vector&,double* der=NULL); - double evaluateGaussian(const vector&, const Gaussian&,double* der=NULL); - double getGaussianNormalization( const Gaussian& ); - vector getGaussianSupport(const Gaussian&); - bool scanOneHill(IFile *ifile, vector &v, vector ¢er, vector &sigma, double &height, bool &multivariate); + double getBias(const std::vector&); + double getBiasAndDerivatives(const std::vector&, std::vector&); + double evaluateGaussian(const std::vector&, const Gaussian&); + double evaluateGaussianAndDerivatives(const std::vector&, const Gaussian&,std::vector&,std::vector&); + double getGaussianNormalization(const Gaussian&); + std::vector getGaussianSupport(const Gaussian&); + bool scanOneHill(IFile* ifile, std::vector& v, std::vector& center, std::vector& sigma, double& height, bool& multivariate); void computeReweightingFactor(); double getTransitionBarrierBias(); - void updateFrequencyAdaptiveStride(); - string fmt; + void updateFrequencyAdaptiveStride(); + void updateNlist(); public: explicit MetaD(const ActionOptions&); @@ -489,11 +517,13 @@ void MetaD::registerKeywords(Keywords& keys) { keys.addOutputComponent("rbias","CALC_RCT","the instantaneous value of the bias normalized using the \\f$c(t)\\f$ reweighting factor [rbias=bias-rct]." "This component can be used to obtain a reweighted histogram."); keys.addOutputComponent("rct","CALC_RCT","the reweighting factor \\f$c(t)\\f$."); - keys.addOutputComponent("work","default","accumulator for work"); + keys.addOutputComponent("work","CALC_WORK","accumulator for work"); keys.addOutputComponent("acc","ACCELERATION","the metadynamics acceleration factor"); keys.addOutputComponent("maxbias", "CALC_MAX_BIAS", "the maximum of the metadynamics V(s, t)"); keys.addOutputComponent("transbias", "CALC_TRANSITION_BIAS", "the metadynamics transition bias V*(t)"); keys.addOutputComponent("pace","FREQUENCY_ADAPTIVE","the hill addition frequency when employing frequency adaptive metadynamics"); + keys.addOutputComponent("nlker","NLIST","number of hills in the neighbor list"); + keys.addOutputComponent("nlsteps","NLIST","number of steps from last neighbor list update"); keys.use("ARG"); keys.add("compulsory","SIGMA","the widths of the Gaussian hills"); keys.add("compulsory","PACE","the frequency for hill addition"); @@ -501,6 +531,7 @@ void MetaD::registerKeywords(Keywords& keys) { keys.add("optional","HEIGHT","the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR are given"); keys.add("optional","FMT","specify format for HILLS files (useful for decrease the number of digits in regtests)"); keys.add("optional","BIASFACTOR","use well tempered metadynamics and use this bias factor. Please note you must also specify temp"); + keys.addFlag("CALC_WORK",false,"calculate the work done by the bias between each update"); keys.add("optional","RECT","list of bias factors for all the replicas"); keys.add("optional","DAMPFACTOR","damp hills with exp(-max(V)/(\\f$k_B\\f$T*DAMPFACTOR)"); for (size_t i = 0; i < n_tempering_options_; i++) { @@ -509,29 +540,31 @@ void MetaD::registerKeywords(Keywords& keys) { keys.add("optional","TARGET","target to a predefined distribution"); keys.add("optional","TEMP","the system temperature - this is only needed if you are doing well-tempered metadynamics"); keys.add("optional","TAU","in well tempered metadynamics, sets height to (\\f$k_B \\Delta T\\f$*pace*timestep)/tau"); - keys.add("optional","GRID_MIN","the lower bounds for the grid"); - keys.add("optional","GRID_MAX","the upper bounds for the grid"); - keys.add("optional","GRID_BIN","the number of bins for the grid"); - keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)"); keys.addFlag("CALC_RCT",false,"calculate the \\f$c(t)\\f$ reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]." "This method is not compatible with metadynamics not on a grid."); keys.add("optional","RCT_USTRIDE","the update stride for calculating the \\f$c(t)\\f$ reweighting factor." "The default 1, so \\f$c(t)\\f$ is updated every time the bias is updated."); + keys.add("optional","GRID_MIN","the lower bounds for the grid"); + keys.add("optional","GRID_MAX","the upper bounds for the grid"); + keys.add("optional","GRID_BIN","the number of bins for the grid"); + keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)"); keys.addFlag("GRID_SPARSE",false,"use a sparse grid to store hills"); keys.addFlag("GRID_NOSPLINE",false,"don't use spline interpolation with grids"); keys.add("optional","GRID_WSTRIDE","write the grid to a file every N steps"); keys.add("optional","GRID_WFILE","the file on which to write the grid"); keys.add("optional","GRID_RFILE","a grid file from which the bias should be read at the initial step of the simulation"); keys.addFlag("STORE_GRIDS",false,"store all the grid files the calculation generates. They will be deleted if this keyword is not present"); + keys.addFlag("NLIST",false,"Use neighbor list for kernels summation, faster but experimental"); + keys.add("optional", "NLIST_PARAMETERS","(default=6.,0.5) the two cutoff parameters for the Gaussians neighbor list"); keys.add("optional","ADAPTIVE","use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or time step dimensions"); + keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds "); + keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds "); keys.add("optional","WALKERS_ID", "walker id"); keys.add("optional","WALKERS_N", "number of walkers"); keys.add("optional","WALKERS_DIR", "shared directory with the hills files from all the walkers"); keys.add("optional","WALKERS_RSTRIDE","stride for reading hills files"); - keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force."); - keys.add("optional","SIGMA_MAX","the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds "); - keys.add("optional","SIGMA_MIN","the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds "); keys.addFlag("WALKERS_MPI",false,"Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR"); + keys.add("optional","INTERVAL","one dimensional lower and upper limits, outside the limits the system will not feel the biasing force."); keys.addFlag("FLYING_GAUSSIAN",false,"Switch on flying Gaussian method, must be used with WALKERS_MPI"); keys.addFlag("ACCELERATION",false,"Set to TRUE if you want to compute the metadynamics acceleration factor."); keys.add("optional","ACCELERATION_RFILE","a data file from which the acceleration should be read at the initial step of the simulation"); @@ -557,40 +590,40 @@ void MetaD::registerTemperingKeywords(const std::string &name_stem, const std::s MetaD::MetaD(const ActionOptions& ao): PLUMED_BIAS_INIT(ao), -// Grid stuff initialization - wgridstride_(0), grid_(false), -// Metadynamics basic parameters - height0_(std::numeric_limits::max()), biasf_(-1.0), dampfactor_(0.0), - tt_specs_(false, "TT", "Transition Tempered", -1.0, 0.0, 1.0), kbt_(0.0), - stride_(0), welltemp_(false), -// frequency adaptive - current_stride(0), - freq_adaptive_(false), - fa_update_frequency_(0), - fa_max_stride_(0), - fa_min_acceleration_(1.0), -// Other stuff + stride_(0), + calc_work_(false), + welltemp_(false), + biasf_(-1.0), + isFirstStep_(true), + height0_(std::numeric_limits::max()), adaptive_(FlexibleBin::none), -// Multiple walkers initialization + grid_(false), + wgridstride_(0), mw_n_(1), mw_dir_(""), mw_id_(0), mw_rstride_(1), - walkers_mpi(false), mpi_nw_(0), mpi_mw_(0), -// Flying Gaussian - flying(false), - acceleration(false), acc(0.0), acc_restart_mean_(0.0), + walkers_mpi_(false), mpi_nw_(0), mpi_mw_(0), + flying_(false), + acceleration_(false), acc_(0.0), acc_restart_mean_(0.0), calc_max_bias_(false), max_bias_(0.0), calc_transition_bias_(false), transition_bias_(0.0), -// Interval initialization + dampfactor_(0.0), + tt_specs_(false, "TT", "Transition Tempered", -1.0, 0.0, 1.0), + current_stride_(0), + freq_adaptive_(false), + fa_update_frequency_(0), + fa_max_stride_(0), + fa_min_acceleration_(1.0), uppI_(-1), lowI_(-1), doInt_(false), - isFirstStep(true), calc_rct_(false), reweight_factor_(0.0), rct_ustride_(1), work_(0), - last_step_warn_grid(0) + nlist_(false), + nlist_update_(false), + nlist_steps_(0) { // parse the flexible hills - string adaptiveoption; + std::string adaptiveoption; adaptiveoption="NONE"; parse("ADAPTIVE",adaptiveoption); if(adaptiveoption=="GEOM") { @@ -605,7 +638,7 @@ MetaD::MetaD(const ActionOptions& ao): error("I do not know this type of adaptive scheme"); } - parse("FMT",fmt); + parse("FMT",fmt_); // parse the sigma parseVector("SIGMA",sigma0_); @@ -640,14 +673,15 @@ MetaD::MetaD(const ActionOptions& ao): for(unsigned i=0; i(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_); + flexbin_=Tools::make_unique(adaptive_,this,sigma0_[0],sigma0min_,sigma0max_); } + // note: HEIGHT is not compulsory, since one could use the TAU keyword, see below parse("HEIGHT",height0_); parse("PACE",stride_); if(stride_<=0 ) error("frequency for hill addition is nonsensical"); - current_stride = stride_; - string hillsfname="HILLS"; + current_stride_ = stride_; + std::string hillsfname="HILLS"; parse("FILE",hillsfname); // Manually set to calculate special bias quantities @@ -681,6 +715,8 @@ MetaD::MetaD(const ActionOptions& ao): if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, with damped metad you must specify it using TEMP"); } + parseFlag("CALC_WORK",calc_work_); + // Set transition tempering parameters. // Transition wells are read later via calc_transition_bias_. readTemperingSpecs(tt_specs_); @@ -689,7 +725,7 @@ MetaD::MetaD(const ActionOptions& ao): // If any previous option specified to calculate a transition bias, // now read the transition wells for that quantity. if (calc_transition_bias_) { - vector tempcoords(getNumberOfArguments()); + std::vector tempcoords(getNumberOfArguments()); for (unsigned i = 0; ; i++) { if (!parseNumberedVector("TRANSITIONWELL", i, tempcoords) ) break; if (tempcoords.size() != getNumberOfArguments()) { @@ -719,21 +755,21 @@ MetaD::MetaD(const ActionOptions& ao): } // Grid Stuff - vector gmin(getNumberOfArguments()); + std::vector gmin(getNumberOfArguments()); parseVector("GRID_MIN",gmin); if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN"); - vector gmax(getNumberOfArguments()); + std::vector gmax(getNumberOfArguments()); parseVector("GRID_MAX",gmax); if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX"); - vector gbin(getNumberOfArguments()); - vector gspacing; + std::vector gbin(getNumberOfArguments()); + std::vector gspacing; parseVector("GRID_BIN",gbin); if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN"); parseVector("GRID_SPACING",gspacing); if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING"); if(gmin.size()!=gmax.size()) error("GRID_MAX and GRID_MIN should be either present or absent"); - if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present"); - if(gbin.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN should be present"); + if(gspacing.size()!=0 && gmin.size()==0) error("If GRID_SPACING is present also GRID_MIN and GRID_MAX should be present"); + if(gbin.size()!=0 && gmin.size()==0) error("If GRID_BIN is present also GRID_MIN and GRID_MAX should be present"); if(gmin.size()!=0) { if(gbin.size()==0 && gspacing.size()==0) { if(adaptive_==FlexibleBin::none) { @@ -743,7 +779,7 @@ MetaD::MetaD(const ActionOptions& ao): for(unsigned i=0; i0) grid_=true; + bool sparsegrid=false; parseFlag("GRID_SPARSE",sparsegrid); bool nospline=false; parseFlag("GRID_NOSPLINE",nospline); bool spline=!nospline; - if(gbin.size()>0) {grid_=true;} parse("GRID_WSTRIDE",wgridstride_); - string gridfilename_; + std::string gridfilename_; parse("GRID_WFILE",gridfilename_); parseFlag("STORE_GRIDS",storeOldGrids_); if(grid_ && gridfilename_.length()>0) { if(wgridstride_==0 ) error("frequency with which to output grid not specified use GRID_WSTRIDE"); } - if(grid_ && wgridstride_>0) { if(gridfilename_.length()==0) error("grid filename not specified use GRID_WFILE"); } - string gridreadfilename_; + + std::string gridreadfilename_; parse("GRID_RFILE",gridreadfilename_); if(!grid_&&gridfilename_.length()> 0) error("To write a grid you need first to define it!"); if(!grid_&&gridreadfilename_.length()>0) error("To read a grid you need first to define it!"); + /*setup neighbor list stuff*/ + parseFlag("NLIST", nlist_); + nlist_center_.resize(getNumberOfArguments()); + nlist_dev2_.resize(getNumberOfArguments()); + if(nlist_&&grid_) error("NLIST and GRID cannot be combined!"); + std::vector nlist_param; + parseVector("NLIST_PARAMETERS",nlist_param); + if(nlist_param.size()==0) + { + nlist_param_[0]=6.0;//*DP2CUTOFF -> max distance of neighbors + nlist_param_[1]=0.5;//*nlist_dev2_[i] -> condition for rebuilding + } + else + { + plumed_massert(nlist_param.size()==2,"two cutoff parameters are needed for the neighbor list"); + plumed_massert(nlist_param[0]>1.0,"the first of NLIST_PARAMETERS must be greater than 1. The smaller the first, the smaller should be the second as well"); + const double min_PARAM_1=(1.-1./std::sqrt(nlist_param[0]/2))+0.16; + plumed_massert(nlist_param[1]>0,"the second of NLIST_PARAMETERS must be greater than 0"); + plumed_massert(nlist_param[1]<=min_PARAM_1,"the second of NLIST_PARAMETERS must be smaller to avoid systematic errors. Largest suggested value is: 1.16-1/sqrt(PARAM_0/2) = "+std::to_string(min_PARAM_1)); + nlist_param_[0]=nlist_param[0]; + nlist_param_[1]=nlist_param[1]; + } + // Reweighting factor rct parseFlag("CALC_RCT",calc_rct_); - if (calc_rct_) - plumed_massert(grid_,"CALC_RCT is supported only if bias is on a grid"); + if (calc_rct_) plumed_massert(grid_,"CALC_RCT is supported only if bias is on a grid"); parse("RCT_USTRIDE",rct_ustride_); if(dampfactor_>0.0) { @@ -804,13 +863,13 @@ MetaD::MetaD(const ActionOptions& ao): parse("WALKERS_RSTRIDE",mw_rstride_); // MPI version - parseFlag("WALKERS_MPI",walkers_mpi); + parseFlag("WALKERS_MPI",walkers_mpi_); // Flying Gaussian - parseFlag("FLYING_GAUSSIAN", flying); + parseFlag("FLYING_GAUSSIAN", flying_); // Inteval keyword - vector tmpI(2); + std::vector tmpI(2); parseVector("INTERVAL",tmpI); if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL"); else if(tmpI.size()==2) { @@ -822,11 +881,10 @@ MetaD::MetaD(const ActionOptions& ao): doInt_=true; } - acceleration=false; - parseFlag("ACCELERATION",acceleration); + parseFlag("ACCELERATION",acceleration_); // Check for a restart acceleration if acceleration is active. - string acc_rfilename; - if (acceleration) { + std::string acc_rfilename; + if(acceleration_) { parse("ACCELERATION_RFILE", acc_rfilename); } @@ -836,7 +894,7 @@ MetaD::MetaD(const ActionOptions& ao): fa_update_frequency_=0; parse("FA_UPDATE_FREQUENCY",fa_update_frequency_); if(fa_update_frequency_!=0 && !freq_adaptive_) { - plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); + plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive METAD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); } if(fa_update_frequency_==0 && freq_adaptive_) { fa_update_frequency_=stride_; @@ -845,13 +903,13 @@ MetaD::MetaD(const ActionOptions& ao): fa_max_stride_=0; parse("FA_MAX_PACE",fa_max_stride_); if(fa_max_stride_!=0 && !freq_adaptive_) { - plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); + plumed_merror("It doesn't make sense to use the FA_MAX_PACE keyword if frequency adaptive METAD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); } // fa_min_acceleration_=1.0; parse("FA_MIN_ACCELERATION",fa_min_acceleration_); if(fa_min_acceleration_!=1.0 && !freq_adaptive_) { - plumed_merror("It doesn't make sense to use the FA_MIN_ACCELERATION keyword if frequency adaptive MetaD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); + plumed_merror("It doesn't make sense to use the FA_MIN_ACCELERATION keyword if frequency adaptive METAD hasn't been activated by using the FREQUENCY_ADAPTIVE flag"); } checkRead(); @@ -868,6 +926,7 @@ MetaD::MetaD(const ActionOptions& ao): log.printf(" Hills relaxation time (tau) %f\n",tau); log.printf(" KbT %f\n",kbt_); } + // Transition tempered metadynamics options if (tt_specs_.is_active) { logTemperingSpecs(tt_specs_); @@ -895,14 +954,14 @@ MetaD::MetaD(const ActionOptions& ao): if (tt_specs_.is_active) n_active++; // Find the greatest alpha. double greatest_alpha = 0.0; - if (welltemp_) greatest_alpha = max(greatest_alpha, 1.0); - if (dampfactor_ > 0.0) greatest_alpha = max(greatest_alpha, 1.0); - if (tt_specs_.is_active) greatest_alpha = max(greatest_alpha, tt_specs_.alpha); + if (welltemp_) greatest_alpha = std::max(greatest_alpha, 1.0); + if (dampfactor_ > 0.0) greatest_alpha = std::max(greatest_alpha, 1.0); + if (tt_specs_.is_active) greatest_alpha = std::max(greatest_alpha, tt_specs_.alpha); // Find the least alpha. double least_alpha = 1.0; - if (welltemp_) least_alpha = min(least_alpha, 1.0); - if (dampfactor_ > 0.0) least_alpha = min(least_alpha, 1.0); - if (tt_specs_.is_active) least_alpha = min(least_alpha, tt_specs_.alpha); + if (welltemp_) least_alpha = std::min(least_alpha, 1.0); + if (dampfactor_ > 0.0) least_alpha = std::min(least_alpha, 1.0); + if (tt_specs_.is_active) least_alpha = std::min(least_alpha, tt_specs_.alpha); // Find the inverse harmonic average of the delta T parameters for all // of the temperings with the greatest alpha values. double total_governing_deltaT_inv = 0.0; @@ -925,6 +984,7 @@ MetaD::MetaD(const ActionOptions& ao): } if(doInt_) log.printf(" Upper and Lower limits boundaries for the bias are activated at %f - %f\n", lowI_, uppI_); + if(grid_) { log.printf(" Grid min"); for(unsigned i=0; i1) { - if(walkers_mpi) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers"); + if(walkers_mpi_) error("MPI version of multiple walkers is not compatible with filesystem version of multiple walkers"); log.printf(" %d multiple walkers active\n",mw_n_); log.printf(" walker id %d\n",mw_id_); log.printf(" reading stride %d\n",mw_rstride_); if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str()); } else { - if(walkers_mpi) { + if(walkers_mpi_) { log.printf(" Multiple walkers active using MPI communnication\n"); if(mw_dir_!="")log.printf(" directory with hills files %s\n",mw_dir_.c_str()); if(comm.Get_rank()==0) { @@ -962,20 +1022,28 @@ MetaD::MetaD(const ActionOptions& ao): } } - if(flying) { - if(!walkers_mpi) error("Flying Gaussian method must be used with MPI version of multiple walkers"); + if(flying_) { + if(!walkers_mpi_) error("Flying Gaussian method must be used with MPI version of multiple walkers"); log.printf(" Flying Gaussian method with %d walkers active\n",mpi_nw_); } + if(nlist_) { + addComponent("nlker"); + componentIsNotPeriodic("nlker"); + addComponent("nlsteps"); + componentIsNotPeriodic("nlsteps"); + } + if(calc_rct_) { addComponent("rbias"); componentIsNotPeriodic("rbias"); addComponent("rct"); componentIsNotPeriodic("rct"); log.printf(" The c(t) reweighting factor will be calculated every %u hills\n",rct_ustride_); getPntrToComponent("rct")->set(reweight_factor_); } - addComponent("work"); componentIsNotPeriodic("work"); - if(acceleration) { + if(calc_work_) { addComponent("work"); componentIsNotPeriodic("work"); } + + if(acceleration_) { if (kbt_ == 0.0) { error("The calculation of the acceleration works only if simulation temperature has been defined"); } @@ -986,7 +1054,8 @@ MetaD::MetaD(const ActionOptions& ao): if (acc_rfilename.length() == 0) { getPntrToComponent("acc")->set(1.0); if(getRestart()) { - log.printf(" WARNING: calculating the acceleration factor in a restarted run without reading in the previous value will most likely lead to incorrect results. You should use the ACCELERATION_RFILE keyword.\n"); + log.printf(" WARNING: calculating the acceleration factor in a restarted run without reading in the previous value will most likely lead to incorrect results.\n"); + log.printf(" You should use the ACCELERATION_RFILE keyword.\n"); } // Otherwise, read and set the restart value. } else { @@ -1017,12 +1086,14 @@ MetaD::MetaD(const ActionOptions& ao): log.printf(" initial acceleration factor read from file %s: value of %f at time %f\n",acc_rfilename.c_str(),acc_rmean,acc_rtime); } } + if (calc_max_bias_) { if (!grid_) error("Calculating the maximum bias on the fly works only with a grid"); log.printf(" calculation on the fly of the maximum bias max(V(s,t)) \n"); addComponent("maxbias"); componentIsNotPeriodic("maxbias"); } + if (calc_transition_bias_) { if (!grid_) error("Calculating the transition bias on the fly works only with a grid"); log.printf(" calculation on the fly of the transition bias V*(t)\n"); @@ -1049,16 +1120,17 @@ MetaD::MetaD(const ActionOptions& ao): } if(freq_adaptive_) { - if(!acceleration) { + if(!acceleration_) { plumed_merror("Frequency adaptive metadynamics only works if the calculation of the acceleration factor is enabled with the ACCELERATION keyword\n"); } - if(walkers_mpi) { + if(walkers_mpi_) { plumed_merror("Combining frequency adaptive metadynamics with MPI multiple walkers is not allowed"); } log.printf(" Frequency adaptive metadynamics enabled\n"); if(getRestart() && acc_rfilename.length() == 0) { - log.printf(" WARNING: using the frequency adaptive scheme in a restarted run without reading in the previous value of the acceleration factor will most likely lead to incorrect results. You should use the ACCELERATION_RFILE keyword.\n"); + log.printf(" WARNING: using the frequency adaptive scheme in a restarted run without reading in the previous value of the acceleration factor will most likely lead to incorrect results.\n"); + log.printf(" You should use the ACCELERATION_RFILE keyword.\n"); } log.printf(" The frequency for hill addition will change dynamically based on the metadynamics acceleration factor\n"); log.printf(" The hill addition frequency will be updated every %d steps\n",fa_update_frequency_); @@ -1072,126 +1144,145 @@ MetaD::MetaD(const ActionOptions& ao): updateFrequencyAdaptiveStride(); } - // for performance - dp_=Tools::make_unique(getNumberOfArguments()); - // initializing and checking grid + bool restartedFromGrid=false; // restart from grid file if(grid_) { - // check for mesh and sigma size - for(unsigned i=0; i0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n"; + if(!(gridreadfilename_.length()>0)) { + // check for mesh and sigma size + for(unsigned i=0; i0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width (SIGMA) can produce artifacts\n"; + } else { + if(sigma0min_[i]<0.) error("When using ADAPTIVE Gaussians on a grid SIGMA_MIN must be specified"); + if(mesh>0.5*sigma0min_[i]) log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians (SIGMA_MIN) \n"; + } + } + std::string funcl=getLabel() + ".bias"; + if(!sparsegrid) {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} + else {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} + std::vector actualmin=BiasGrid_->getMin(); + std::vector actualmax=BiasGrid_->getMax(); + for(unsigned i=0; i0.5*sigma0min_[i]||sigma0min_[i]<0.) log<<" WARNING: to use a METAD with a GRID and ADAPTIVE you need to set a Grid Spacing larger than half of the Gaussians \n"; + error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!"); } - } - std::string funcl=getLabel() + ".bias"; - if(!sparsegrid) {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} - else {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} - std::vector actualmin=BiasGrid_->getMin(); - std::vector actualmax=BiasGrid_->getMax(); - for(unsigned i=0; igetDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments"); + for(unsigned i=0; iisPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias"); + double a, b; + Tools::convert(gmin[i],a); + Tools::convert(gmax[i],b); + double mesh=(b-a)/((double)gbin[i]); + if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n"; + } + log.printf(" Restarting from %s\n",gridreadfilename_.c_str()); + if(getRestart()) restartedFromGrid=true; } } - // restart from external grid - bool restartedFromGrid=false; - if(gridreadfilename_.length()>0) { - // read the grid in input, find the keys - IFile gridfile; - gridfile.link(*this); - if(gridfile.FileExist(gridreadfilename_)) { - gridfile.open(gridreadfilename_); - } else { - error("The GRID file you want to read: " + gridreadfilename_ + ", cannot be found!"); - } - std::string funcl=getLabel() + ".bias"; - BiasGrid_=GridBase::create(funcl, getArguments(), gridfile, gmin, gmax, gbin, sparsegrid, spline, true); - if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments"); - for(unsigned i=0; iisPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias"); - double a, b; - Tools::convert(gmin[i],a); - Tools::convert(gmax[i],b); - double mesh=(b-a)/((double)gbin[i]); - if(mesh>0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n"; - } - log.printf(" Restarting from %s:",gridreadfilename_.c_str()); - if(getRestart()) restartedFromGrid=true; + // if we are restarting from GRID and using WALKERS_MPI we can check that all walkers have actually read the grid + if(getRestart()&&walkers_mpi_) { + std::vector restarted(mpi_nw_,0); + if(comm.Get_rank()==0) multi_sim_comm.Allgather(int(restartedFromGrid), restarted); + comm.Bcast(restarted,0); + int result = std::accumulate(restarted.begin(),restarted.end(),0); + if(result!=0&&result!=mpi_nw_) error("in this WALKERS_MPI run some replica have restarted from GRID while other do not!"); } - // initializing and checking grid - if(grid_&&!(gridreadfilename_.length()>0)) { - // check for adaptive and sigma_min - if(sigma0min_.size()==0&&adaptive_!=FlexibleBin::none) error("When using Adaptive Gaussians on a grid SIGMA_MIN must be specified"); - // check for mesh and sigma size - for(unsigned i=0; i0.5*sigma0_[i]) log<<" WARNING: Using a METAD with a Grid Spacing larger than half of the Gaussians width can produce artifacts\n"; - } - std::string funcl=getLabel() + ".bias"; - if(!sparsegrid) {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} - else {BiasGrid_=Tools::make_unique(funcl,getArguments(),gmin,gmax,gbin,spline,true);} - std::vector actualmin=BiasGrid_->getMin(); - std::vector actualmax=BiasGrid_->getMax(); - for(unsigned i=0; i1) { - stringstream out; out << i; + std::stringstream out; out << i; fname = mw_dir_+"/"+hillsfname+"."+out.str(); - } else if(walkers_mpi) { + } else if(walkers_mpi_) { fname = mw_dir_+"/"+hillsfname; } else { fname = hillsfname; } } else { if(mw_n_>1) { - stringstream out; out << i; + std::stringstream out; out << i; fname = hillsfname+"."+out.str(); } else { fname = hillsfname; } } - ifiles.emplace_back(Tools::make_unique()); + ifiles_.emplace_back(Tools::make_unique()); // this is just a shortcut pointer to the last element: - IFile *ifile = ifiles.back().get(); - ifilesnames.push_back(fname); + IFile *ifile = ifiles_.back().get(); + ifilesnames_.push_back(fname); ifile->link(*this); if(ifile->FileExist(fname)) { ifile->open(fname); if(getRestart()&&!restartedFromGrid) { - log.printf(" Restarting from %s:",ifilesnames[i].c_str()); - readGaussians(ifiles[i].get()); + log.printf(" Restarting from %s:",ifilesnames_[i].c_str()); + readGaussians(ifiles_[i].get()); + restartedFromHills=true; } - ifiles[i]->reset(false); + ifiles_[i]->reset(false); // close only the walker own hills file for later writing - if(i==mw_id_) ifiles[i]->close(); + if(i==mw_id_) ifiles_[i]->close(); } else { // in case a file does not exist and we are restarting, complain that the file was not found - if(getRestart()) log<<" WARNING: restart file "< restarted(mpi_nw_,0); + if(comm.Get_rank()==0) multi_sim_comm.Allgather(int(restartedFromHills), restarted); + comm.Bcast(restarted,0); + int result = std::accumulate(restarted.begin(),restarted.end(),0); + if(result!=0&&result!=mpi_nw_) error("in this WALKERS_MPI run some replica have restarted from FILE while other do not!"); + } + comm.Barrier(); // this barrier is needed when using walkers_mpi // to be sure that all files have been read before @@ -1200,7 +1291,8 @@ MetaD::MetaD(const ActionOptions& ao): // it would introduce troubles when using replicas without METAD // (e.g. in bias exchange with a neutral replica) // see issue #168 on github - if(comm.Get_rank()==0 && walkers_mpi) multi_sim_comm.Barrier(); + if(comm.Get_rank()==0 && walkers_mpi_) multi_sim_comm.Barrier(); + if(targetfilename_.length()>0) { IFile gridfile; gridfile.open(targetfilename_); std::string funcl=getLabel() + ".target"; @@ -1211,22 +1303,26 @@ MetaD::MetaD(const ActionOptions& ao): } } - // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills - if(getRestart() && calc_rct_) computeReweightingFactor(); - // Calculate all special bias quantities desired if restarting with nonzero bias. - if(getRestart() && calc_max_bias_) { - max_bias_ = BiasGrid_->getMaxValue(); - getPntrToComponent("maxbias")->set(max_bias_); - } - if(getRestart() && calc_transition_bias_) { - transition_bias_ = getTransitionBarrierBias(); - getPntrToComponent("transbias")->set(transition_bias_); + if(getRestart()) { + // if this is a restart the neighbor list should be immediately updated + if(nlist_) nlist_update_=true; + // Calculate the Tiwary-Parrinello reweighting factor if we are restarting from previous hills + if(calc_rct_) computeReweightingFactor(); + // Calculate all special bias quantities desired if restarting with nonzero bias. + if(calc_max_bias_) { + max_bias_ = BiasGrid_->getMaxValue(); + getPntrToComponent("maxbias")->set(max_bias_); + } + if(calc_transition_bias_) { + transition_bias_ = getTransitionBarrierBias(); + getPntrToComponent("transbias")->set(transition_bias_); + } } // open grid file for writing if(wgridstride_>0) { gridfile_.link(*this); - if(walkers_mpi) { + if(walkers_mpi_) { int r=0; if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank(); comm.Bcast(r,0); @@ -1239,16 +1335,16 @@ MetaD::MetaD(const ActionOptions& ao): // open hills file for writing hillsOfile_.link(*this); - if(walkers_mpi) { + if(walkers_mpi_) { int r=0; if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank(); comm.Bcast(r,0); - if(r>0) ifilesnames[mw_id_]="/dev/null"; + if(r>0) ifilesnames_[mw_id_]="/dev/null"; hillsOfile_.enforceSuffix(""); } if(mw_n_>1) hillsOfile_.enforceSuffix(""); - hillsOfile_.open(ifilesnames[mw_id_]); - if(fmt.length()>0) hillsOfile_.fmtField(fmt); + hillsOfile_.open(ifilesnames_[mw_id_]); + if(fmt_.length()>0) hillsOfile_.fmtField(fmt_); hillsOfile_.addConstantField("multivariate"); hillsOfile_.addConstantField("kerneltype"); if(doInt_) { @@ -1264,7 +1360,7 @@ MetaD::MetaD(const ActionOptions& ao): for(const auto & p : actionSet) if(dynamic_cast(p.get())) { concurrent=true; break; } if(concurrent) log<<" You are using concurrent metadynamics\n"; if(rect_biasf_.size()>0) { - if(walkers_mpi) { + if(walkers_mpi_) { log<<" You are using RECT in its 'altruistic' implementation\n"; }{ log<<" You are using RECT\n"; @@ -1272,38 +1368,29 @@ MetaD::MetaD(const ActionOptions& ao): } log<<" Bibliography "<1||walkers_mpi) log<0) log<0 && walkers_mpi) log<1||walkers_mpi_) log<0) log<0 && walkers_mpi_) log<0) { log< center(ncv); - vector sigma(ncv); + std::vector center(ncv); + std::vector sigma(ncv); double height; int nhills=0; bool multivariate=false; @@ -1344,12 +1432,12 @@ void MetaD::readGaussians(IFile *ifile) std::vector tmpvalues; for(unsigned j=0; jgetName(), false ) ); - while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) { - ; + while(scanOneHill(ifile,tmpvalues,center,sigma,height,multivariate)) + { nhills++; -// note that for gamma=1 we store directly -F - if(welltemp_ && biasf_>1.0) {height*=(biasf_-1.0)/biasf_;} - addGaussian(Gaussian(center,sigma,height,multivariate)); + // note that for gamma=1 we store directly -F + if(welltemp_ && biasf_>1.0) height*=(biasf_-1.0)/biasf_; + addGaussian(Gaussian(multivariate,height,center,sigma)); } log.printf(" %d Gaussians read\n",nhills); } @@ -1398,7 +1486,7 @@ void MetaD::writeGaussian(const Gaussian& hill, OFile&file) file.printField("sigma_"+getPntrToArgument(i)->getName(),hill.sigma[i]); } double height=hill.height; -// note that for gamma=1 we store directly -F + // note that for gamma=1 we store directly -F if(welltemp_ && biasf_>1.0) height*=biasf_/(biasf_-1.0); file.printField("height",height).printField("biasf",biasf_); if(mw_n_>1) file.printField("clock",int(std::time(0))); @@ -1407,46 +1495,52 @@ void MetaD::writeGaussian(const Gaussian& hill, OFile&file) void MetaD::addGaussian(const Gaussian& hill) { - if(!grid_) hills_.push_back(hill); - else { + if(grid_) { size_t ncv=getNumberOfArguments(); - vector nneighb=getGaussianSupport(hill); - vector neighbors=BiasGrid_->getNeighbors(hill.center,nneighb); - vector der(ncv); - vector xx(ncv); + std::vector nneighb=getGaussianSupport(hill); + std::vector neighbors=BiasGrid_->getNeighbors(hill.center,nneighb); + std::vector der(ncv); + std::vector xx(ncv); if(comm.Get_size()==1) { + // for performance reasons and thread safety + std::vector dp(ncv); for(size_t i=0; igetPoint(ineigh,xx); - double bias=evaluateGaussian(xx,hill,&der[0]); + double bias=evaluateGaussianAndDerivatives(xx,hill,der,dp); BiasGrid_->addValueAndDerivatives(ineigh,bias,der); } } else { unsigned stride=comm.Get_size(); unsigned rank=comm.Get_rank(); - vector allder(ncv*neighbors.size(),0.0); - vector allbias(neighbors.size(),0.0); + std::vector allder(ncv*neighbors.size(),0.0); + std::vector n_der(ncv,0.0); + std::vector allbias(neighbors.size(),0.0); + // for performance reasons and thread safety + std::vector dp(ncv); for(unsigned i=rank; igetPoint(ineigh,xx); - allbias[i]=evaluateGaussian(xx,hill,&allder[ncv*i]); + allbias[i]=evaluateGaussianAndDerivatives(xx,hill,n_der,dp); + for(unsigned j=0; jaddValueAndDerivatives(ineigh,allbias[i],der); } } - } + } else hills_.push_back(hill); } -vector MetaD::getGaussianSupport(const Gaussian& hill) +std::vector MetaD::getGaussianSupport(const Gaussian& hill) { - vector nneigh; - vector cutoff; + std::vector nneigh; + std::vector cutoff; unsigned ncv=getNumberOfArguments(); // traditional or flexible hill? @@ -1464,7 +1558,7 @@ vector MetaD::getGaussianSupport(const Gaussian& hill) Matrix myinv(ncv,ncv); Invert(mymatrix,myinv); Matrix myautovec(ncv,ncv); - vector myautoval(ncv); //should I take this or their square root? + std::vector myautoval(ncv); //should I take this or their square root? diagMat(myinv,myautoval,myautovec); double maxautoval=0.; unsigned ind_maxautoval; ind_maxautoval=ncv; @@ -1472,11 +1566,11 @@ vector MetaD::getGaussianSupport(const Gaussian& hill) if(myautoval[i]>maxautoval) {maxautoval=myautoval[i]; ind_maxautoval=i;} } for(unsigned i=0; i MetaD::getGaussianSupport(const Gaussian& hill) return nneigh; } -double MetaD::getBiasAndDerivatives(const vector& cv, double* der) +double MetaD::getBias(const std::vector& cv) { double bias=0.0; - if(!grid_) { - if(hills_.size()>10000 && (getStep()-last_step_warn_grid)>10000) { - std::string msg; - Tools::convert(hills_.size(),msg); - msg="You have accumulated "+msg+" hills, you should enable GRIDs to avoid serious performance hits"; - warning(msg); - last_step_warn_grid=getStep(); - } + if(grid_) bias = BiasGrid_->getValue(cv); + else { + unsigned nt=OpenMP::getNumThreads(); unsigned stride=comm.Get_size(); unsigned rank=comm.Get_rank(); - for(unsigned i=rank; i& cv, std::vector& der) +{ + unsigned ncv=getNumberOfArguments(); + double bias=0.0; + if(grid_) { + std::vector vder(ncv); + bias=BiasGrid_->getValueAndDerivatives(cv,vder); + for(unsigned i=0; i vder(getNumberOfArguments()); - bias=BiasGrid_->getValueAndDerivatives(cv,vder); - for(unsigned i=0; i dp(ncv); + for(unsigned i=rank; i omp_deriv(ncv,0.); + // for performance reasons and thread safety + std::vector dp(ncv); + #pragma omp for reduction(+:bias) nowait + for(unsigned i=rank; igetValue(cv); + if(hills_.size()<2*nt*stride||nt==1) { + // for performance reasons and thread safety + std::vector dp(ncv); + for(unsigned i=rank; i omp_deriv(ncv,0.); + // for performance reasons and thread safety + std::vector dp(ncv); + #pragma omp for reduction(+:bias) nowait + for(unsigned i=rank; i(ncv)/2.0); + return norm*std::pow(2*pi,static_cast(ncv)/2.0); } -double MetaD::evaluateGaussian(const vector& cv, const Gaussian& hill, double* der) +double MetaD::evaluateGaussian(const std::vector& cv, const Gaussian& hill) { - double dp2=0.0; - double bias=0.0; + unsigned ncv=cv.size(); + // I use a pointer here because cv is const (and should be const) // but when using doInt it is easier to locally replace cv[0] with // the upper/lower limit in case it is out of range + double tmpcv[1]; const double *pcv=NULL; // pointer to cv - double tmpcv[1]; // tmp array with cv (to be used with doInt_) - if(cv.size()>0) pcv=&cv[0]; + if(ncv>0) pcv=&cv[0]; if(doInt_) { - plumed_assert(cv.size()==1); + plumed_assert(ncv==1); tmpcv[0]=cv[0]; if(cv[0]uppI_) tmpcv[0]=uppI_; pcv=&(tmpcv[0]); } + + double dp2=0.0; if(hill.multivariate) { unsigned k=0; - unsigned ncv=cv.size(); // recompose the full sigma from the upper diag cholesky Matrix mymatrix(ncv,ncv); for(unsigned i=0; i& cv, const Gaussian& hill, d k++; } } - for(unsigned i=0; i& cv, const Gaussian& hill, d } } } + } else { + for(unsigned i=0; i& cv, const Gaussian& hill, std::vector& der, std::vector& dp_) +{ + unsigned ncv=cv.size(); + + // I use a pointer here because cv is const (and should be const) + // but when using doInt it is easier to locally replace cv[0] with + // the upper/lower limit in case it is out of range + const double *pcv=NULL; // pointer to cv + double tmpcv[1]; // tmp array with cv (to be used with doInt_) + if(ncv>0) pcv=&cv[0]; + if(doInt_) { + plumed_assert(ncv==1); + tmpcv[0]=cv[0]; + if(cv[0]uppI_) tmpcv[0]=uppI_; + pcv=&(tmpcv[0]); + } + + bool int_der=false; + if(doInt_) { + if(cv[0]uppI_) int_der=true; + } + + double dp2=0.0; + double bias=0.0; + if(hill.multivariate) { + unsigned k=0; + // recompose the full sigma from the upper diag cholesky + Matrix mymatrix(ncv,ncv); + for(unsigned i=0; iuppI_) for(unsigned i=0; i& cv) +double MetaD::getHeight(const std::vector& cv) { double height=height0_; if(welltemp_) { - double vbias = getBiasAndDerivatives(cv); + double vbias = getBias(cv); if(biasf_>1.0) { - height = height0_*exp(-vbias/(kbt_*(biasf_-1.0))); + height = height0_*std::exp(-vbias/(kbt_*(biasf_-1.0))); } else { // notice that if gamma=1 we store directly -F - height = height0_*exp(-vbias/kbt_); + height = height0_*std::exp(-vbias/kbt_); } } if(dampfactor_>0.0) { plumed_assert(BiasGrid_); double m=BiasGrid_->getMaxValue(); - height*=exp(-m/(kbt_*(dampfactor_))); + height*=std::exp(-m/(kbt_*(dampfactor_))); } if (tt_specs_.is_active) { double vbarrier = transition_bias_; @@ -1649,16 +1860,17 @@ double MetaD::getHeight(const vector& cv) } if(TargetGrid_) { double f=TargetGrid_->getValue(cv)-TargetGrid_->getMaxValue(); - height*=exp(f/kbt_); + height*=std::exp(f/kbt_); } return height; } -void MetaD::temperHeight(double &height, const TemperingSpecs &t_specs, const double tempering_bias) { +void MetaD::temperHeight(double& height, const TemperingSpecs& t_specs, const double tempering_bias) +{ if (t_specs.alpha == 1.0) { - height *= exp(-max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0))); + height *= std::exp(-std::max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0))); } else { - height *= pow(1 + (1 - t_specs.alpha) / t_specs.alpha * max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)), - t_specs.alpha / (1 - t_specs.alpha)); + height *= std::pow(1 + (1 - t_specs.alpha) / t_specs.alpha * std::max(0.0, tempering_bias - t_specs.threshold) / (kbt_ * (t_specs.biasf - 1.0)), - t_specs.alpha / (1 - t_specs.alpha)); } } @@ -1669,78 +1881,81 @@ void MetaD::calculate() if(adaptive_==FlexibleBin::diffusion && getExchangeStep()) error("ADAPTIVE=DIFF is not compatible with replica exchange"); const unsigned ncv=getNumberOfArguments(); - vector cv(ncv); - auto der = Tools::make_unique(ncv); - for(unsigned i=0; i cv(ncv); + for(unsigned i=0; inlist_param_[1]) {nlist_update_=true; break;} + } + } + if(nlist_update_) updateNlist(); } + double ene = 0.; + std::vector der(ncv,0.); + if(biasf_!=1.0) ene = getBiasAndDerivatives(cv,der); setBias(ene); + for(unsigned i=0; iset(work_); if(calc_rct_) getPntrToComponent("rbias")->set(ene - reweight_factor_); // calculate the acceleration factor - if(acceleration&&!isFirstStep) { - acc += static_cast(getStride()) * exp(ene/(kbt_)); - const double mean_acc = acc/((double) getStep()); + if(acceleration_&&!isFirstStep_) { + acc_ += static_cast(getStride()) * std::exp(ene/(kbt_)); + const double mean_acc = acc_/((double) getStep()); getPntrToComponent("acc")->set(mean_acc); - } else if (acceleration && isFirstStep && acc_restart_mean_ > 0.0) { - acc = acc_restart_mean_ * static_cast(getStep()); + } else if (acceleration_ && isFirstStep_ && acc_restart_mean_ > 0.0) { + acc_ = acc_restart_mean_ * static_cast(getStep()); if(freq_adaptive_) { // has to be done here if restarting, as the acc is not defined before updateFrequencyAdaptiveStride(); } } - - getPntrToComponent("work")->set(work_); - // set Forces - for(unsigned i=0; i cv(getNumberOfArguments()); - vector thissigma; - bool multivariate; - +void MetaD::update() +{ // adding hills criteria (could be more complex though) bool nowAddAHill; - if(getStep()%current_stride==0 && !isFirstStep )nowAddAHill=true; + if(getStep()%current_stride_==0 && !isFirstStep_) nowAddAHill=true; else { nowAddAHill=false; - isFirstStep=false; + isFirstStep_=false; } - for(unsigned i=0; i cv(ncv); + for(unsigned i=0; iupdate(nowAddAHill); + flexbin_->update(nowAddAHill); multivariate=true; - } else { - multivariate=false; } + std::vector thissigma; if(nowAddAHill) { // add a Gaussian double height=getHeight(cv); // returns upper diagonal inverse - if(adaptive_!=FlexibleBin::none) thissigma=flexbin->getInverseMatrix(); + if(adaptive_!=FlexibleBin::none) thissigma=flexbin_->getInverseMatrix(); // returns normal sigma else thissigma=sigma0_; // In case we use walkers_mpi, it is now necessary to communicate with other replicas. - if(walkers_mpi) { + if(walkers_mpi_) { // Allocate arrays to store all walkers hills - std::vector all_cv(mpi_nw_*cv.size(),0.0); + std::vector all_cv(mpi_nw_*ncv,0.0); std::vector all_sigma(mpi_nw_*thissigma.size(),0.0); std::vector all_height(mpi_nw_,0.0); std::vector all_multivariate(mpi_nw_,0); @@ -1748,7 +1963,7 @@ void MetaD::update() { // Communicate (only root) multi_sim_comm.Allgather(cv,all_cv); multi_sim_comm.Allgather(thissigma,all_sigma); -// notice that if gamma=1 we store directly -F so this scaling is not necessary: + // notice that if gamma=1 we store directly -F so this scaling is not necessary: multi_sim_comm.Allgather(height*(biasf_>1.0?biasf_/(biasf_-1.0):1.0),all_height); multi_sim_comm.Allgather(int(multivariate),all_multivariate); } @@ -1759,43 +1974,42 @@ void MetaD::update() { comm.Bcast(all_multivariate,0); // Flying Gaussian - if (flying) { + if (flying_) { hills_.clear(); comm.Barrier(); } for(unsigned i=0; i cv_now(cv.size()); + std::vector cv_now(ncv); std::vector sigma_now(thissigma.size()); - for(unsigned j=0; j1.0?(biasf_-1.0)/biasf_:1.0),all_multivariate[i]); + // notice that if gamma=1 we store directly -F so this scaling is not necessary: + double fact=(biasf_>1.0?(biasf_-1.0)/biasf_:1.0); + Gaussian newhill=Gaussian(all_multivariate[i],all_height[i]*fact,cv_now,sigma_now); addGaussian(newhill); - - // Flying Gaussian - if (!flying) { - writeGaussian(newhill,hillsOfile_); - } - + if(!flying_) writeGaussian(newhill,hillsOfile_); } } else { - Gaussian newhill=Gaussian(cv,thissigma,height,multivariate); + Gaussian newhill=Gaussian(multivariate,height,cv,thissigma); addGaussian(newhill); - // print on HILLS file writeGaussian(newhill,hillsOfile_); } - } -// this should be outside of the if block in case -// mw_rstride_ is not a multiple of stride_ - if(mw_n_>1 && getStep()%mw_rstride_==0) { - hillsOfile_.flush(); + // this is to update the hills neighbor list + if(nlist_) nlist_update_=true; } - double vbias1=getBiasAndDerivatives(cv); - work_+=vbias1-vbias; + // this should be outside of the if block in case + // mw_rstride_ is not a multiple of stride_ + if(mw_n_>1 && getStep()%mw_rstride_==0) hillsOfile_.flush(); + + if(calc_work_) { + if(nlist_) updateNlist(); + double vbias1=getBias(cv); + work_+=vbias1-vbias; + } // dump grid on file if(wgridstride_>0&&(getStep()%wgridstride_==0||getCPT())) { @@ -1806,7 +2020,7 @@ void MetaD::update() { // this will overwrite previously written grids else { int r = 0; - if(walkers_mpi) { + if(walkers_mpi_) { if(comm.Get_rank()==0) r=multi_sim_comm.Get_rank(); comm.Bcast(r,0); } @@ -1828,20 +2042,23 @@ void MetaD::update() { // don't read your own Gaussians if(i==mw_id_) continue; // if the file is not open yet - if(!(ifiles[i]->isOpen())) { + if(!(ifiles_[i]->isOpen())) { // check if it exists now and open it! - if(ifiles[i]->FileExist(ifilesnames[i])) { - ifiles[i]->open(ifilesnames[i]); - ifiles[i]->reset(false); + if(ifiles_[i]->FileExist(ifilesnames_[i])) { + ifiles_[i]->open(ifilesnames_[i]); + ifiles_[i]->reset(false); } // otherwise read the new Gaussians } else { - log.printf(" Reading hills from %s:",ifilesnames[i].c_str()); - readGaussians(ifiles[i].get()); - ifiles[i]->reset(false); + log.printf(" Reading hills from %s:",ifilesnames_[i].c_str()); + readGaussians(ifiles_[i].get()); + ifiles_[i]->reset(false); } } + // this is to update the hills neighbor list + if(nlist_) nlist_update_=true; } + // Recalculate special bias quantities whenever the bias has been changed by the update. bool bias_has_changed = (nowAddAHill || (mw_n_ > 1 && getStep() % mw_rstride_ == 0)); if (calc_rct_ && bias_has_changed && getStep()%(stride_*rct_ustride_)==0) computeReweightingFactor(); @@ -1858,16 +2075,15 @@ void MetaD::update() { if(freq_adaptive_ && getStep()%fa_update_frequency_==0) { updateFrequencyAdaptiveStride(); } - } -/// takes a pointer to the file and a template string with values v and gives back the next center, sigma and height -bool MetaD::scanOneHill(IFile *ifile, vector &tmpvalues, vector ¢er, vector &sigma, double &height, bool &multivariate) +/// takes a pointer to the file and a template std::string with values v and gives back the next center, sigma and height +bool MetaD::scanOneHill(IFile* ifile, std::vector& tmpvalues, std::vector& center, std::vector& sigma, double& height, bool& multivariate) { double dummy; multivariate=false; if(ifile->scanField("time",dummy)) { - unsigned ncv; ncv=tmpvalues.size(); + unsigned ncv=tmpvalues.size(); for(unsigned i=0; iscanField( &tmpvalues[i] ); if( tmpvalues[i].isPeriodic() && ! getPntrToArgument(i)->isPeriodic() ) { @@ -1955,20 +2171,18 @@ void MetaD::computeReweightingFactor() Z_0+=std::exp(minusBetaF*(val-max_bias_)); Z_V+=std::exp(minusBetaFplusV*(val-max_bias_)); } - if (stride>1) { - comm.Sum(Z_0); - comm.Sum(Z_V); - } + comm.Sum(Z_0); + comm.Sum(Z_V); reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_; getPntrToComponent("rct")->set(reweight_factor_); } -double MetaD::getTransitionBarrierBias() { - +double MetaD::getTransitionBarrierBias() +{ // If there is only one well of interest, return the bias at that well point. if (transitionwells_.size() == 1) { - double tb_bias = getBiasAndDerivatives(transitionwells_[0], NULL); + double tb_bias = getBias(transitionwells_[0]); return tb_bias; // Otherwise, check for the least barrier bias between all pairs of wells. @@ -1983,8 +2197,8 @@ double MetaD::getTransitionBarrierBias() { // transitionwell_[1] is sampled. } else { double least_transition_bias; - vector sink = transitionwells_[0]; - vector source = transitionwells_[1]; + std::vector sink = transitionwells_[0]; + std::vector source = transitionwells_[1]; least_transition_bias = BiasGrid_->findMaximalPathMinimum(source, sink); for (unsigned i = 2; i < transitionwells_.size(); i++) { if (least_transition_bias == 0.0) { @@ -1998,27 +2212,81 @@ double MetaD::getTransitionBarrierBias() { } } -void MetaD::updateFrequencyAdaptiveStride() { +void MetaD::updateFrequencyAdaptiveStride() +{ plumed_massert(freq_adaptive_,"should only be used if frequency adaptive metadynamics is enabled"); - plumed_massert(acceleration,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated"); - const double mean_acc = acc/((double) getStep()); + plumed_massert(acceleration_,"frequency adaptive metadynamics can only be used if the acceleration factor is calculated"); + const double mean_acc = acc_/((double) getStep()); int tmp_stride= stride_*floor((mean_acc/fa_min_acceleration_)+0.5); if(mean_acc >= fa_min_acceleration_) { - if(tmp_stride > current_stride) {current_stride = tmp_stride;} + if(tmp_stride > current_stride_) {current_stride_ = tmp_stride;} } - if(fa_max_stride_!=0 && current_stride>fa_max_stride_) { - current_stride=fa_max_stride_; + if(fa_max_stride_!=0 && current_stride_>fa_max_stride_) { + current_stride_=fa_max_stride_; } - getPntrToComponent("pace")->set(current_stride); + getPntrToComponent("pace")->set(current_stride_); } bool MetaD::checkNeedsGradients()const { if(adaptive_==FlexibleBin::geometry) { - if(getStep()%stride_==0 && !isFirstStep) return true; + if(getStep()%stride_==0 && !isFirstStep_) return true; else return false; } else return false; } +void MetaD::updateNlist() +{ + // no need to check for neighbors + if(hills_.size()==0) return; + + // here we generate the neighbor list + nlist_hills_.clear(); + std::vector local_flat_nl; + unsigned nt=OpenMP::getNumThreads(); + if(hills_.size()<2*nt) nt=1; + #pragma omp parallel num_threads(nt) + { + std::vector private_flat_nl; + #pragma omp for nowait + for(unsigned k=0; k dev2; + dev2.resize(getNumberOfArguments(),0); + for(unsigned k=0; k0.) nlist_dev2_[i]=dev2[i]/static_cast(nlist_hills_.size()); + else nlist_dev2_[i]=hills_.back().sigma[i]*hills_.back().sigma[i]; + } + + // we are done + getPntrToComponent("nlker")->set(nlist_hills_.size()); + getPntrToComponent("nlsteps")->set(nlist_steps_); + nlist_steps_=0; + nlist_update_=false; +} + } }