From 3024e138fd3e77105f0e92fd9ae4f857b07c4ad8 Mon Sep 17 00:00:00 2001 From: Glen Hocky Date: Thu, 14 May 2020 15:10:15 -0400 Subject: [PATCH] New module (FISST) for Infinite Switch Simulated Tempering in Force (#570) * Add code and documentation for FISST; to install, use ./configure --enable-modules=+eds * fixed error message for when no temperature set * added regtests for fisst module, covering a standard case, and an example of a different initial weights distribution * add fisst copyright file, current arxiv reference, and added citations to docs * ran header.sh script on fisst module * made extra library in FISST code use PLMD namespace, to resolve cppcheck error * fix error in include command causing plumedcheck to freeze --- CHANGES/v2.7.md | 4 +- regtest/fisst/Makefile | 2 + regtest/fisst/rt-fisst-basic/Makefile | 1 + regtest/fisst/rt-fisst-basic/config | 5 + .../rt-fisst-basic/observable.out.reference | 11 + regtest/fisst/rt-fisst-basic/plumed.dat | 5 + .../rt-fisst-basic/restart.out.reference | 11 + regtest/fisst/rt-fisst-weightdist/Makefile | 1 + regtest/fisst/rt-fisst-weightdist/config | 5 + .../observable.out.reference | 11 + regtest/fisst/rt-fisst-weightdist/plumed.dat | 5 + .../rt-fisst-weightdist/restart.out.reference | 11 + src/fisst/COPYRIGHT | 14 + src/fisst/FISST.cpp | 661 ++++++++++++++++++ src/fisst/Makefile | 3 + src/fisst/legendre_rule_fast.cpp | 565 +++++++++++++++ src/fisst/legendre_rule_fast.h | 44 ++ src/fisst/module.type | 1 + user-doc/FISSTMOD.md | 26 + user-doc/bibliography.bib | 59 +- 20 files changed, 1416 insertions(+), 29 deletions(-) create mode 100644 regtest/fisst/Makefile create mode 100644 regtest/fisst/rt-fisst-basic/Makefile create mode 100644 regtest/fisst/rt-fisst-basic/config create mode 100644 regtest/fisst/rt-fisst-basic/observable.out.reference create mode 100644 regtest/fisst/rt-fisst-basic/plumed.dat create mode 100644 regtest/fisst/rt-fisst-basic/restart.out.reference create mode 100644 regtest/fisst/rt-fisst-weightdist/Makefile create mode 100644 regtest/fisst/rt-fisst-weightdist/config create mode 100644 regtest/fisst/rt-fisst-weightdist/observable.out.reference create mode 100644 regtest/fisst/rt-fisst-weightdist/plumed.dat create mode 100644 regtest/fisst/rt-fisst-weightdist/restart.out.reference create mode 100644 src/fisst/COPYRIGHT create mode 100644 src/fisst/FISST.cpp create mode 100644 src/fisst/Makefile create mode 100644 src/fisst/legendre_rule_fast.cpp create mode 100644 src/fisst/legendre_rule_fast.h create mode 100644 src/fisst/module.type create mode 100644 user-doc/FISSTMOD.md diff --git a/CHANGES/v2.7.md b/CHANGES/v2.7.md index ef6c2b8b78..f9f52f2f91 100644 --- a/CHANGES/v2.7.md +++ b/CHANGES/v2.7.md @@ -8,9 +8,9 @@ Changes from version 2.6 which are relevant for users: - A new Funnel module by Stefano Raniolo and Vittorio Limongelli - \ref FUNNEL_PS - \ref FUNNEL + - A new Infinite Switch Simulated Tempering in Force module by Glen Hocky + - \ref FISST For developers: - small fix in `Plumed.h` too avoid unique global symbols (see \issue{549}) - - diff --git a/regtest/fisst/Makefile b/regtest/fisst/Makefile new file mode 100644 index 0000000000..42480767ae --- /dev/null +++ b/regtest/fisst/Makefile @@ -0,0 +1,2 @@ +include ../scripts/module.make + diff --git a/regtest/fisst/rt-fisst-basic/Makefile b/regtest/fisst/rt-fisst-basic/Makefile new file mode 100644 index 0000000000..3703b27cea --- /dev/null +++ b/regtest/fisst/rt-fisst-basic/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/fisst/rt-fisst-basic/config b/regtest/fisst/rt-fisst-basic/config new file mode 100644 index 0000000000..f6211702f8 --- /dev/null +++ b/regtest/fisst/rt-fisst-basic/config @@ -0,0 +1,5 @@ +type=driver +plumed_modules=fisst +arg="--plumed plumed.dat --ixyz trajectory.xyz --kt 0.593" +extra_files="../../trajectories/trajectory.xyz" + diff --git a/regtest/fisst/rt-fisst-basic/observable.out.reference b/regtest/fisst/rt-fisst-basic/observable.out.reference new file mode 100644 index 0000000000..4568c9f492 --- /dev/null +++ b/regtest/fisst/rt-fisst-basic/observable.out.reference @@ -0,0 +1,11 @@ +#! FIELDS time nsamples d1 d1_fbar d1_f0 d1_ow0 d1_f1 d1_ow1 d1_f2 d1_ow2 d1_f3 d1_ow3 d1_f4 d1_ow4 d1_f5 d1_ow5 d1_f6 d1_ow6 d1_f7 d1_ow7 d1_f8 d1_ow8 d1_f9 d1_ow9 d1_f10 d1_ow10 d1_f11 d1_ow11 d1_f12 d1_ow12 d1_f13 d1_ow13 d1_f14 d1_ow14 d1_f15 d1_ow15 d1_f16 d1_ow16 d1_f17 d1_ow17 d1_f18 d1_ow18 d1_f19 d1_ow19 d1_f20 d1_ow20 d1_f21 d1_ow21 d1_f22 d1_ow22 d1_f23 d1_ow23 d1_f24 d1_ow24 d1_f25 d1_ow25 d1_f26 d1_ow26 d1_f27 d1_ow27 d1_f28 d1_ow28 d1_f29 d1_ow29 d1_f30 d1_ow30 +#! SET kbt 1.42e-01 +#! SET n_interpolation 31 +#! SET period 1 +#! SET min_force 0.00e+00 +#! SET max_force 1.44e+00 + 0.00e+00 0 1.26e+01 0.00e+00 2.10e-03 1.00e+00 1.10e-02 1.00e+00 2.70e-02 1.00e+00 4.99e-02 1.00e+00 7.94e-02 1.00e+00 1.15e-01 1.00e+00 1.57e-01 1.00e+00 2.05e-01 1.00e+00 2.57e-01 1.00e+00 3.14e-01 1.00e+00 3.76e-01 1.00e+00 4.40e-01 1.00e+00 5.08e-01 1.00e+00 5.77e-01 1.00e+00 6.48e-01 1.00e+00 7.20e-01 1.00e+00 7.92e-01 1.00e+00 8.63e-01 1.00e+00 9.32e-01 1.00e+00 1.00e+00 1.00e+00 1.06e+00 1.00e+00 1.13e+00 1.00e+00 1.18e+00 1.00e+00 1.24e+00 1.00e+00 1.28e+00 1.00e+00 1.32e+00 1.00e+00 1.36e+00 1.00e+00 1.39e+00 1.00e+00 1.41e+00 1.00e+00 1.43e+00 1.00e+00 1.44e+00 1.00e+00 + 1.00e+03 2 1.32e+01 7.20e-01 2.10e-03 1.31e+02 1.10e-02 1.31e+02 2.70e-02 1.31e+02 4.99e-02 1.31e+02 7.94e-02 1.31e+02 1.15e-01 1.31e+02 1.57e-01 1.31e+02 2.05e-01 1.31e+02 2.57e-01 1.31e+02 3.14e-01 1.31e+02 3.76e-01 1.31e+02 4.40e-01 1.31e+02 5.08e-01 1.31e+02 5.77e-01 1.31e+02 6.48e-01 1.31e+02 7.20e-01 1.31e+02 7.92e-01 1.31e+02 8.63e-01 1.31e+02 9.32e-01 1.31e+02 1.00e+00 1.31e+02 1.06e+00 1.31e+02 1.13e+00 1.31e+02 1.18e+00 1.31e+02 1.24e+00 1.31e+02 1.28e+00 1.31e+02 1.32e+00 1.31e+02 1.36e+00 1.31e+02 1.39e+00 1.31e+02 1.41e+00 1.31e+02 1.43e+00 1.31e+02 1.44e+00 1.31e+02 + 2.00e+03 3 1.39e+01 8.67e-01 2.10e-03 3.15e+01 1.10e-02 3.28e+01 2.70e-02 3.53e+01 4.99e-02 3.90e+01 7.94e-02 4.44e+01 1.15e-01 5.17e+01 1.57e-01 6.12e+01 2.05e-01 7.32e+01 2.57e-01 8.77e+01 3.14e-01 1.05e+02 3.76e-01 1.23e+02 4.40e-01 1.42e+02 5.08e-01 1.60e+02 5.77e-01 1.77e+02 6.48e-01 1.90e+02 7.20e-01 2.00e+02 7.92e-01 2.08e+02 8.63e-01 2.14e+02 9.32e-01 2.18e+02 1.00e+00 2.21e+02 1.06e+00 2.22e+02 1.13e+00 2.24e+02 1.18e+00 2.25e+02 1.24e+00 2.25e+02 1.28e+00 2.26e+02 1.32e+00 2.26e+02 1.36e+00 2.26e+02 1.39e+00 2.26e+02 1.41e+00 2.26e+02 1.43e+00 2.27e+02 1.44e+00 2.27e+02 + 3.00e+03 4 1.48e+01 9.49e-01 2.10e-03 3.81e+00 1.10e-02 4.17e+00 2.70e-02 4.91e+00 4.99e-02 6.19e+00 7.94e-02 8.30e+00 1.15e-01 1.18e+01 1.57e-01 1.74e+01 2.05e-01 2.67e+01 2.57e-01 4.14e+01 3.14e-01 6.38e+01 3.76e-01 9.53e+01 4.40e-01 1.35e+02 5.08e-01 1.79e+02 5.77e-01 2.21e+02 6.48e-01 2.57e+02 7.20e-01 2.85e+02 7.92e-01 3.05e+02 8.63e-01 3.19e+02 9.32e-01 3.29e+02 1.00e+00 3.35e+02 1.06e+00 3.39e+02 1.13e+00 3.42e+02 1.18e+00 3.44e+02 1.24e+00 3.45e+02 1.28e+00 3.46e+02 1.32e+00 3.46e+02 1.36e+00 3.47e+02 1.39e+00 3.47e+02 1.41e+00 3.47e+02 1.43e+00 3.47e+02 1.44e+00 3.47e+02 + 4.00e+03 5 1.49e+01 7.93e-01 2.10e-03 4.69e+01 1.10e-02 5.09e+01 2.70e-02 5.87e+01 4.99e-02 7.10e+01 7.94e-02 8.88e+01 1.15e-01 1.12e+02 1.57e-01 1.40e+02 2.05e-01 1.69e+02 2.57e-01 1.95e+02 3.14e-01 2.14e+02 3.76e-01 2.28e+02 4.40e-01 2.36e+02 5.08e-01 2.41e+02 5.77e-01 2.45e+02 6.48e-01 2.47e+02 7.20e-01 2.48e+02 7.92e-01 2.49e+02 8.63e-01 2.49e+02 9.32e-01 2.50e+02 1.00e+00 2.50e+02 1.06e+00 2.51e+02 1.13e+00 2.51e+02 1.18e+00 2.51e+02 1.24e+00 2.51e+02 1.28e+00 2.52e+02 1.32e+00 2.52e+02 1.36e+00 2.52e+02 1.39e+00 2.52e+02 1.41e+00 2.52e+02 1.43e+00 2.52e+02 1.44e+00 2.52e+02 diff --git a/regtest/fisst/rt-fisst-basic/plumed.dat b/regtest/fisst/rt-fisst-basic/plumed.dat new file mode 100644 index 0000000000..70af2d2243 --- /dev/null +++ b/regtest/fisst/rt-fisst-basic/plumed.dat @@ -0,0 +1,5 @@ +UNITS LENGTH=A TIME=fs ENERGY=kcal/mol +d1: DISTANCE ATOMS=1,2 + +#0 pN to 100 pN +f: FISST ARG=d1 PERIOD=1 MIN_FORCE=0 MAX_FORCE=1.44 NINTERPOLATE=31 OUT_RESTART=restart.out OUT_OBSERVABLE=observable.out RESTART_FMT=%4.2e diff --git a/regtest/fisst/rt-fisst-basic/restart.out.reference b/regtest/fisst/rt-fisst-basic/restart.out.reference new file mode 100644 index 0000000000..59c3c1c2f4 --- /dev/null +++ b/regtest/fisst/rt-fisst-basic/restart.out.reference @@ -0,0 +1,11 @@ +#! FIELDS time nsamples d1 d1_f0 d1_g0 d1_w0 d1_z0 d1_f1 d1_g1 d1_w1 d1_z1 d1_f2 d1_g2 d1_w2 d1_z2 d1_f3 d1_g3 d1_w3 d1_z3 d1_f4 d1_g4 d1_w4 d1_z4 d1_f5 d1_g5 d1_w5 d1_z5 d1_f6 d1_g6 d1_w6 d1_z6 d1_f7 d1_g7 d1_w7 d1_z7 d1_f8 d1_g8 d1_w8 d1_z8 d1_f9 d1_g9 d1_w9 d1_z9 d1_f10 d1_g10 d1_w10 d1_z10 d1_f11 d1_g11 d1_w11 d1_z11 d1_f12 d1_g12 d1_w12 d1_z12 d1_f13 d1_g13 d1_w13 d1_z13 d1_f14 d1_g14 d1_w14 d1_z14 d1_f15 d1_g15 d1_w15 d1_z15 d1_f16 d1_g16 d1_w16 d1_z16 d1_f17 d1_g17 d1_w17 d1_z17 d1_f18 d1_g18 d1_w18 d1_z18 d1_f19 d1_g19 d1_w19 d1_z19 d1_f20 d1_g20 d1_w20 d1_z20 d1_f21 d1_g21 d1_w21 d1_z21 d1_f22 d1_g22 d1_w22 d1_z22 d1_f23 d1_g23 d1_w23 d1_z23 d1_f24 d1_g24 d1_w24 d1_z24 d1_f25 d1_g25 d1_w25 d1_z25 d1_f26 d1_g26 d1_w26 d1_z26 d1_f27 d1_g27 d1_w27 d1_z27 d1_f28 d1_g28 d1_w28 d1_z28 d1_f29 d1_g29 d1_w29 d1_z29 d1_f30 d1_g30 d1_w30 d1_z30 +#! SET kbt 1.42e-01 +#! SET n_interpolation 31 +#! SET period 1 +#! SET min_force 0.00e+00 +#! SET max_force 1.44e+00 + 0.00e+00 0 1.26e+01 2.10e-03 5.38e-03 6.94e-01 1.44e+00 1.10e-02 1.25e-02 6.94e-01 1.44e+00 2.70e-02 1.94e-02 6.94e-01 1.44e+00 4.99e-02 2.62e-02 6.94e-01 1.44e+00 7.94e-02 3.28e-02 6.94e-01 1.44e+00 1.15e-01 3.90e-02 6.94e-01 1.44e+00 1.57e-01 4.48e-02 6.94e-01 1.44e+00 2.05e-01 5.01e-02 6.94e-01 1.44e+00 2.57e-01 5.50e-02 6.94e-01 1.44e+00 3.14e-01 5.93e-02 6.94e-01 1.44e+00 3.76e-01 6.31e-02 6.94e-01 1.44e+00 4.40e-01 6.62e-02 6.94e-01 1.44e+00 5.08e-01 6.86e-02 6.94e-01 1.44e+00 5.77e-01 7.04e-02 6.94e-01 1.44e+00 6.48e-01 7.14e-02 6.94e-01 1.44e+00 7.20e-01 7.18e-02 6.94e-01 1.44e+00 7.92e-01 7.14e-02 6.94e-01 1.44e+00 8.63e-01 7.04e-02 6.94e-01 1.44e+00 9.32e-01 6.86e-02 6.94e-01 1.44e+00 1.00e+00 6.62e-02 6.94e-01 1.44e+00 1.06e+00 6.31e-02 6.94e-01 1.44e+00 1.13e+00 5.93e-02 6.94e-01 1.44e+00 1.18e+00 5.50e-02 6.94e-01 1.44e+00 1.24e+00 5.01e-02 6.94e-01 1.44e+00 1.28e+00 4.48e-02 6.94e-01 1.44e+00 1.32e+00 3.90e-02 6.94e-01 1.44e+00 1.36e+00 3.28e-02 6.94e-01 1.44e+00 1.39e+00 2.62e-02 6.94e-01 1.44e+00 1.41e+00 1.94e-02 6.94e-01 1.44e+00 1.43e+00 1.25e-02 6.94e-01 1.44e+00 1.44e+00 5.38e-03 6.94e-01 1.44e+00 + 1.00e+03 2 1.32e+01 2.10e-03 5.38e-03 7.65e+01 9.95e-05 1.10e-02 1.25e-02 3.34e+01 2.28e-04 2.70e-02 1.94e-02 7.56e+00 1.01e-03 4.99e-02 2.62e-02 9.03e-01 8.43e-03 7.94e-02 3.28e-02 5.80e-02 1.31e-01 1.15e-01 3.90e-02 2.07e-03 3.69e+00 1.57e-01 4.48e-02 4.20e-05 1.81e+02 2.05e-01 5.01e-02 5.08e-07 1.50e+04 2.57e-01 5.50e-02 3.82e-09 1.99e+06 3.14e-01 5.93e-02 1.87e-11 4.06e+08 3.76e-01 6.31e-02 6.31e-14 1.21e+11 4.40e-01 6.62e-02 1.55e-16 4.92e+13 5.08e-01 6.86e-02 2.93e-19 2.60e+16 5.77e-01 7.04e-02 4.55e-22 1.67e+19 6.48e-01 7.14e-02 6.21e-25 1.23e+22 7.20e-01 7.18e-02 7.93e-28 9.60e+24 7.92e-01 7.14e-02 1.01e-30 7.52e+27 8.63e-01 7.04e-02 1.38e-33 5.51e+30 9.32e-01 6.86e-02 2.15e-36 3.54e+33 1.00e+00 6.62e-02 4.07e-39 1.87e+36 1.06e+00 6.31e-02 9.96e-42 7.64e+38 1.13e+00 5.93e-02 3.36e-44 2.27e+41 1.18e+00 5.50e-02 1.64e-46 4.63e+43 1.24e+00 5.01e-02 1.24e-48 6.16e+45 1.28e+00 4.48e-02 1.50e-50 5.09e+47 1.32e+00 3.90e-02 3.04e-52 2.50e+49 1.36e+00 3.28e-02 1.08e-53 7.03e+50 1.39e+00 2.62e-02 6.96e-55 1.09e+52 1.41e+00 1.94e-02 8.32e-56 9.15e+52 1.43e+00 1.25e-02 1.88e-56 4.04e+53 1.44e+00 5.38e-03 8.21e-57 9.26e+53 + 2.00e+03 3 1.39e+01 2.10e-03 5.38e-03 7.70e+01 7.70e-05 1.10e-02 1.25e-02 3.34e+01 1.78e-04 2.70e-02 1.94e-02 7.46e+00 7.95e-04 4.99e-02 2.62e-02 8.73e-01 6.79e-03 7.94e-02 3.28e-02 5.46e-02 1.09e-01 1.15e-01 3.90e-02 1.86e-03 3.18e+00 1.57e-01 4.48e-02 3.59e-05 1.65e+02 2.05e-01 5.01e-02 4.03e-07 1.47e+04 2.57e-01 5.50e-02 2.74e-09 2.16e+06 3.14e-01 5.93e-02 1.18e-11 5.03e+08 3.76e-01 6.31e-02 3.37e-14 1.76e+11 4.40e-01 6.62e-02 6.75e-17 8.78e+13 5.08e-01 6.86e-02 1.00e-19 5.90e+16 5.77e-01 7.04e-02 1.19e-22 5.00e+19 6.48e-01 7.14e-02 1.19e-25 4.98e+22 7.20e-01 7.18e-02 1.09e-28 5.42e+25 7.92e-01 7.14e-02 9.88e-32 6.00e+28 8.63e-01 7.04e-02 9.47e-35 6.26e+31 9.32e-01 6.86e-02 1.03e-37 5.73e+34 1.00e+00 6.62e-02 1.38e-40 4.29e+37 1.06e+00 6.31e-02 2.42e-43 2.45e+40 1.13e+00 5.93e-02 5.91e-46 1.00e+43 1.18e+00 5.50e-02 2.14e-48 2.77e+45 1.24e+00 5.01e-02 1.22e-50 4.87e+47 1.28e+00 4.48e-02 1.14e-52 5.18e+49 1.32e+00 3.90e-02 1.86e-54 3.18e+51 1.36e+00 3.28e-02 5.48e-56 1.08e+53 1.39e+00 2.62e-02 3.01e-57 1.97e+54 1.41e+00 1.94e-02 3.18e-58 1.86e+55 1.43e+00 1.25e-02 6.62e-59 8.95e+55 1.44e+00 5.38e-03 2.75e-59 2.15e+56 + 3.00e+03 4 1.48e+01 2.10e-03 5.38e-03 7.71e+01 5.84e-05 1.10e-02 1.25e-02 3.34e+01 1.35e-04 2.70e-02 1.94e-02 7.44e+00 6.05e-04 4.99e-02 2.62e-02 8.68e-01 5.18e-03 7.94e-02 3.28e-02 5.39e-02 8.35e-02 1.15e-01 3.90e-02 1.82e-03 2.47e+00 1.57e-01 4.48e-02 3.45e-05 1.30e+02 2.05e-01 5.01e-02 3.76e-07 1.20e+04 2.57e-01 5.50e-02 2.44e-09 1.84e+06 3.14e-01 5.93e-02 9.75e-12 4.62e+08 3.76e-01 6.31e-02 2.48e-14 1.82e+11 4.40e-01 6.62e-02 4.19e-17 1.08e+14 5.08e-01 6.86e-02 4.95e-20 9.09e+16 5.77e-01 7.04e-02 4.39e-23 1.03e+20 6.48e-01 7.14e-02 3.16e-26 1.43e+23 7.20e-01 7.18e-02 2.01e-29 2.24e+26 7.92e-01 7.14e-02 1.24e-32 3.64e+29 8.63e-01 7.04e-02 8.00e-36 5.62e+32 9.32e-01 6.86e-02 5.91e-39 7.62e+35 1.00e+00 6.62e-02 5.38e-42 8.37e+38 1.06e+00 6.31e-02 6.49e-45 6.93e+41 1.13e+00 5.93e-02 1.11e-47 4.04e+44 1.18e+00 5.50e-02 2.90e-50 1.55e+47 1.24e+00 5.01e-02 1.22e-52 3.69e+49 1.28e+00 4.48e-02 8.71e-55 5.17e+51 1.32e+00 3.90e-02 1.11e-56 4.04e+53 1.36e+00 3.28e-02 2.66e-58 1.69e+55 1.39e+00 2.62e-02 1.23e-59 3.66e+56 1.41e+00 1.94e-02 1.14e-60 3.95e+57 1.43e+00 1.25e-02 2.16e-61 2.08e+58 1.44e+00 5.38e-03 8.54e-62 5.27e+58 + 4.00e+03 5 1.49e+01 2.10e-03 5.38e-03 7.86e+01 5.72e-05 1.10e-02 1.25e-02 3.34e+01 1.35e-04 2.70e-02 1.94e-02 7.16e+00 6.28e-04 4.99e-02 2.62e-02 7.83e-01 5.75e-03 7.94e-02 3.28e-02 4.39e-02 1.03e-01 1.15e-01 3.90e-02 1.27e-03 3.53e+00 1.57e-01 4.48e-02 1.94e-05 2.32e+02 2.05e-01 5.01e-02 1.58e-07 2.84e+04 2.57e-01 5.50e-02 7.21e-10 6.25e+06 3.14e-01 5.93e-02 1.93e-12 2.33e+09 3.76e-01 6.31e-02 3.28e-15 1.37e+12 4.40e-01 6.62e-02 3.78e-18 1.19e+15 5.08e-01 6.86e-02 3.21e-21 1.40e+18 5.77e-01 7.04e-02 2.16e-24 2.08e+21 6.48e-01 7.14e-02 1.25e-27 3.60e+24 7.20e-01 7.18e-02 6.67e-31 6.74e+27 7.92e-01 7.14e-02 3.56e-34 1.26e+31 8.63e-01 7.04e-02 2.05e-37 2.20e+34 9.32e-01 6.86e-02 1.36e-40 3.30e+37 1.00e+00 6.62e-02 1.13e-43 3.97e+40 1.06e+00 6.31e-02 1.26e-46 3.56e+43 1.13e+00 5.93e-02 2.02e-49 2.23e+46 1.18e+00 5.50e-02 4.92e-52 9.15e+48 1.24e+00 5.01e-02 1.94e-54 2.31e+51 1.28e+00 4.48e-02 1.32e-56 3.42e+53 1.32e+00 3.90e-02 1.61e-58 2.80e+55 1.36e+00 3.28e-02 3.69e-60 1.22e+57 1.39e+00 2.62e-02 1.66e-61 2.72e+58 1.41e+00 1.94e-02 1.50e-62 3.01e+59 1.43e+00 1.25e-02 2.79e-63 1.61e+60 1.44e+00 5.38e-03 1.09e-63 4.13e+60 diff --git a/regtest/fisst/rt-fisst-weightdist/Makefile b/regtest/fisst/rt-fisst-weightdist/Makefile new file mode 100644 index 0000000000..3703b27cea --- /dev/null +++ b/regtest/fisst/rt-fisst-weightdist/Makefile @@ -0,0 +1 @@ +include ../../scripts/test.make diff --git a/regtest/fisst/rt-fisst-weightdist/config b/regtest/fisst/rt-fisst-weightdist/config new file mode 100644 index 0000000000..f6211702f8 --- /dev/null +++ b/regtest/fisst/rt-fisst-weightdist/config @@ -0,0 +1,5 @@ +type=driver +plumed_modules=fisst +arg="--plumed plumed.dat --ixyz trajectory.xyz --kt 0.593" +extra_files="../../trajectories/trajectory.xyz" + diff --git a/regtest/fisst/rt-fisst-weightdist/observable.out.reference b/regtest/fisst/rt-fisst-weightdist/observable.out.reference new file mode 100644 index 0000000000..4568c9f492 --- /dev/null +++ b/regtest/fisst/rt-fisst-weightdist/observable.out.reference @@ -0,0 +1,11 @@ +#! FIELDS time nsamples d1 d1_fbar d1_f0 d1_ow0 d1_f1 d1_ow1 d1_f2 d1_ow2 d1_f3 d1_ow3 d1_f4 d1_ow4 d1_f5 d1_ow5 d1_f6 d1_ow6 d1_f7 d1_ow7 d1_f8 d1_ow8 d1_f9 d1_ow9 d1_f10 d1_ow10 d1_f11 d1_ow11 d1_f12 d1_ow12 d1_f13 d1_ow13 d1_f14 d1_ow14 d1_f15 d1_ow15 d1_f16 d1_ow16 d1_f17 d1_ow17 d1_f18 d1_ow18 d1_f19 d1_ow19 d1_f20 d1_ow20 d1_f21 d1_ow21 d1_f22 d1_ow22 d1_f23 d1_ow23 d1_f24 d1_ow24 d1_f25 d1_ow25 d1_f26 d1_ow26 d1_f27 d1_ow27 d1_f28 d1_ow28 d1_f29 d1_ow29 d1_f30 d1_ow30 +#! SET kbt 1.42e-01 +#! SET n_interpolation 31 +#! SET period 1 +#! SET min_force 0.00e+00 +#! SET max_force 1.44e+00 + 0.00e+00 0 1.26e+01 0.00e+00 2.10e-03 1.00e+00 1.10e-02 1.00e+00 2.70e-02 1.00e+00 4.99e-02 1.00e+00 7.94e-02 1.00e+00 1.15e-01 1.00e+00 1.57e-01 1.00e+00 2.05e-01 1.00e+00 2.57e-01 1.00e+00 3.14e-01 1.00e+00 3.76e-01 1.00e+00 4.40e-01 1.00e+00 5.08e-01 1.00e+00 5.77e-01 1.00e+00 6.48e-01 1.00e+00 7.20e-01 1.00e+00 7.92e-01 1.00e+00 8.63e-01 1.00e+00 9.32e-01 1.00e+00 1.00e+00 1.00e+00 1.06e+00 1.00e+00 1.13e+00 1.00e+00 1.18e+00 1.00e+00 1.24e+00 1.00e+00 1.28e+00 1.00e+00 1.32e+00 1.00e+00 1.36e+00 1.00e+00 1.39e+00 1.00e+00 1.41e+00 1.00e+00 1.43e+00 1.00e+00 1.44e+00 1.00e+00 + 1.00e+03 2 1.32e+01 7.20e-01 2.10e-03 1.31e+02 1.10e-02 1.31e+02 2.70e-02 1.31e+02 4.99e-02 1.31e+02 7.94e-02 1.31e+02 1.15e-01 1.31e+02 1.57e-01 1.31e+02 2.05e-01 1.31e+02 2.57e-01 1.31e+02 3.14e-01 1.31e+02 3.76e-01 1.31e+02 4.40e-01 1.31e+02 5.08e-01 1.31e+02 5.77e-01 1.31e+02 6.48e-01 1.31e+02 7.20e-01 1.31e+02 7.92e-01 1.31e+02 8.63e-01 1.31e+02 9.32e-01 1.31e+02 1.00e+00 1.31e+02 1.06e+00 1.31e+02 1.13e+00 1.31e+02 1.18e+00 1.31e+02 1.24e+00 1.31e+02 1.28e+00 1.31e+02 1.32e+00 1.31e+02 1.36e+00 1.31e+02 1.39e+00 1.31e+02 1.41e+00 1.31e+02 1.43e+00 1.31e+02 1.44e+00 1.31e+02 + 2.00e+03 3 1.39e+01 8.67e-01 2.10e-03 3.15e+01 1.10e-02 3.28e+01 2.70e-02 3.53e+01 4.99e-02 3.90e+01 7.94e-02 4.44e+01 1.15e-01 5.17e+01 1.57e-01 6.12e+01 2.05e-01 7.32e+01 2.57e-01 8.77e+01 3.14e-01 1.05e+02 3.76e-01 1.23e+02 4.40e-01 1.42e+02 5.08e-01 1.60e+02 5.77e-01 1.77e+02 6.48e-01 1.90e+02 7.20e-01 2.00e+02 7.92e-01 2.08e+02 8.63e-01 2.14e+02 9.32e-01 2.18e+02 1.00e+00 2.21e+02 1.06e+00 2.22e+02 1.13e+00 2.24e+02 1.18e+00 2.25e+02 1.24e+00 2.25e+02 1.28e+00 2.26e+02 1.32e+00 2.26e+02 1.36e+00 2.26e+02 1.39e+00 2.26e+02 1.41e+00 2.26e+02 1.43e+00 2.27e+02 1.44e+00 2.27e+02 + 3.00e+03 4 1.48e+01 9.49e-01 2.10e-03 3.81e+00 1.10e-02 4.17e+00 2.70e-02 4.91e+00 4.99e-02 6.19e+00 7.94e-02 8.30e+00 1.15e-01 1.18e+01 1.57e-01 1.74e+01 2.05e-01 2.67e+01 2.57e-01 4.14e+01 3.14e-01 6.38e+01 3.76e-01 9.53e+01 4.40e-01 1.35e+02 5.08e-01 1.79e+02 5.77e-01 2.21e+02 6.48e-01 2.57e+02 7.20e-01 2.85e+02 7.92e-01 3.05e+02 8.63e-01 3.19e+02 9.32e-01 3.29e+02 1.00e+00 3.35e+02 1.06e+00 3.39e+02 1.13e+00 3.42e+02 1.18e+00 3.44e+02 1.24e+00 3.45e+02 1.28e+00 3.46e+02 1.32e+00 3.46e+02 1.36e+00 3.47e+02 1.39e+00 3.47e+02 1.41e+00 3.47e+02 1.43e+00 3.47e+02 1.44e+00 3.47e+02 + 4.00e+03 5 1.49e+01 7.93e-01 2.10e-03 4.69e+01 1.10e-02 5.09e+01 2.70e-02 5.87e+01 4.99e-02 7.10e+01 7.94e-02 8.88e+01 1.15e-01 1.12e+02 1.57e-01 1.40e+02 2.05e-01 1.69e+02 2.57e-01 1.95e+02 3.14e-01 2.14e+02 3.76e-01 2.28e+02 4.40e-01 2.36e+02 5.08e-01 2.41e+02 5.77e-01 2.45e+02 6.48e-01 2.47e+02 7.20e-01 2.48e+02 7.92e-01 2.49e+02 8.63e-01 2.49e+02 9.32e-01 2.50e+02 1.00e+00 2.50e+02 1.06e+00 2.51e+02 1.13e+00 2.51e+02 1.18e+00 2.51e+02 1.24e+00 2.51e+02 1.28e+00 2.52e+02 1.32e+00 2.52e+02 1.36e+00 2.52e+02 1.39e+00 2.52e+02 1.41e+00 2.52e+02 1.43e+00 2.52e+02 1.44e+00 2.52e+02 diff --git a/regtest/fisst/rt-fisst-weightdist/plumed.dat b/regtest/fisst/rt-fisst-weightdist/plumed.dat new file mode 100644 index 0000000000..8a6017b457 --- /dev/null +++ b/regtest/fisst/rt-fisst-weightdist/plumed.dat @@ -0,0 +1,5 @@ +UNITS LENGTH=A TIME=fs ENERGY=kcal/mol +d1: DISTANCE ATOMS=1,2 + +#0 pN to 100 pN +f: FISST ARG=d1 PERIOD=1 MIN_FORCE=0 MAX_FORCE=1.44 INITIAL_WEIGHT_DIST=EXP INITIAL_WEIGHT_RATE=0.7 NINTERPOLATE=31 OUT_RESTART=restart.out OUT_OBSERVABLE=observable.out RESTART_FMT=%4.2e diff --git a/regtest/fisst/rt-fisst-weightdist/restart.out.reference b/regtest/fisst/rt-fisst-weightdist/restart.out.reference new file mode 100644 index 0000000000..3e998a9cd8 --- /dev/null +++ b/regtest/fisst/rt-fisst-weightdist/restart.out.reference @@ -0,0 +1,11 @@ +#! FIELDS time nsamples d1 d1_f0 d1_g0 d1_w0 d1_z0 d1_f1 d1_g1 d1_w1 d1_z1 d1_f2 d1_g2 d1_w2 d1_z2 d1_f3 d1_g3 d1_w3 d1_z3 d1_f4 d1_g4 d1_w4 d1_z4 d1_f5 d1_g5 d1_w5 d1_z5 d1_f6 d1_g6 d1_w6 d1_z6 d1_f7 d1_g7 d1_w7 d1_z7 d1_f8 d1_g8 d1_w8 d1_z8 d1_f9 d1_g9 d1_w9 d1_z9 d1_f10 d1_g10 d1_w10 d1_z10 d1_f11 d1_g11 d1_w11 d1_z11 d1_f12 d1_g12 d1_w12 d1_z12 d1_f13 d1_g13 d1_w13 d1_z13 d1_f14 d1_g14 d1_w14 d1_z14 d1_f15 d1_g15 d1_w15 d1_z15 d1_f16 d1_g16 d1_w16 d1_z16 d1_f17 d1_g17 d1_w17 d1_z17 d1_f18 d1_g18 d1_w18 d1_z18 d1_f19 d1_g19 d1_w19 d1_z19 d1_f20 d1_g20 d1_w20 d1_z20 d1_f21 d1_g21 d1_w21 d1_z21 d1_f22 d1_g22 d1_w22 d1_z22 d1_f23 d1_g23 d1_w23 d1_z23 d1_f24 d1_g24 d1_w24 d1_z24 d1_f25 d1_g25 d1_w25 d1_z25 d1_f26 d1_g26 d1_w26 d1_z26 d1_f27 d1_g27 d1_w27 d1_z27 d1_f28 d1_g28 d1_w28 d1_z28 d1_f29 d1_g29 d1_w29 d1_z29 d1_f30 d1_g30 d1_w30 d1_z30 +#! SET kbt 1.42e-01 +#! SET n_interpolation 31 +#! SET period 1 +#! SET min_force 0.00e+00 +#! SET max_force 1.44e+00 + 0.00e+00 0 1.26e+01 2.10e-03 5.38e-03 1.10e+00 9.09e-01 1.10e-02 1.25e-02 1.09e+00 9.14e-01 2.70e-02 1.94e-02 1.08e+00 9.25e-01 4.99e-02 2.62e-02 1.06e+00 9.39e-01 7.94e-02 3.28e-02 1.04e+00 9.59e-01 1.15e-01 3.90e-02 1.02e+00 9.83e-01 1.57e-01 4.48e-02 9.87e-01 1.01e+00 2.05e-01 5.01e-02 9.55e-01 1.05e+00 2.57e-01 5.50e-02 9.21e-01 1.09e+00 3.14e-01 5.93e-02 8.84e-01 1.13e+00 3.76e-01 6.31e-02 8.47e-01 1.18e+00 4.40e-01 6.62e-02 8.10e-01 1.23e+00 5.08e-01 6.86e-02 7.73e-01 1.29e+00 5.77e-01 7.04e-02 7.36e-01 1.36e+00 6.48e-01 7.14e-02 7.00e-01 1.43e+00 7.20e-01 7.18e-02 6.66e-01 1.50e+00 7.92e-01 7.14e-02 6.33e-01 1.58e+00 8.63e-01 7.04e-02 6.03e-01 1.66e+00 9.32e-01 6.86e-02 5.74e-01 1.74e+00 1.00e+00 6.62e-02 5.48e-01 1.83e+00 1.06e+00 6.31e-02 5.23e-01 1.91e+00 1.13e+00 5.93e-02 5.01e-01 1.99e+00 1.18e+00 5.50e-02 4.82e-01 2.08e+00 1.24e+00 5.01e-02 4.64e-01 2.15e+00 1.28e+00 4.48e-02 4.49e-01 2.23e+00 1.32e+00 3.90e-02 4.36e-01 2.29e+00 1.36e+00 3.28e-02 4.25e-01 2.35e+00 1.39e+00 2.62e-02 4.17e-01 2.40e+00 1.41e+00 1.94e-02 4.10e-01 2.44e+00 1.43e+00 1.25e-02 4.05e-01 2.47e+00 1.44e+00 5.38e-03 4.03e-01 2.48e+00 + 1.00e+03 2 1.32e+01 2.10e-03 5.38e-03 7.65e+01 9.95e-05 1.10e-02 1.25e-02 3.34e+01 2.28e-04 2.70e-02 1.94e-02 7.56e+00 1.01e-03 4.99e-02 2.62e-02 9.03e-01 8.43e-03 7.94e-02 3.28e-02 5.80e-02 1.31e-01 1.15e-01 3.90e-02 2.07e-03 3.69e+00 1.57e-01 4.48e-02 4.20e-05 1.81e+02 2.05e-01 5.01e-02 5.08e-07 1.50e+04 2.57e-01 5.50e-02 3.82e-09 1.99e+06 3.14e-01 5.93e-02 1.87e-11 4.06e+08 3.76e-01 6.31e-02 6.31e-14 1.21e+11 4.40e-01 6.62e-02 1.55e-16 4.92e+13 5.08e-01 6.86e-02 2.93e-19 2.60e+16 5.77e-01 7.04e-02 4.55e-22 1.67e+19 6.48e-01 7.14e-02 6.21e-25 1.23e+22 7.20e-01 7.18e-02 7.93e-28 9.60e+24 7.92e-01 7.14e-02 1.01e-30 7.52e+27 8.63e-01 7.04e-02 1.38e-33 5.51e+30 9.32e-01 6.86e-02 2.15e-36 3.54e+33 1.00e+00 6.62e-02 4.07e-39 1.87e+36 1.06e+00 6.31e-02 9.96e-42 7.64e+38 1.13e+00 5.93e-02 3.36e-44 2.27e+41 1.18e+00 5.50e-02 1.64e-46 4.63e+43 1.24e+00 5.01e-02 1.24e-48 6.16e+45 1.28e+00 4.48e-02 1.50e-50 5.09e+47 1.32e+00 3.90e-02 3.04e-52 2.50e+49 1.36e+00 3.28e-02 1.08e-53 7.03e+50 1.39e+00 2.62e-02 6.96e-55 1.09e+52 1.41e+00 1.94e-02 8.32e-56 9.15e+52 1.43e+00 1.25e-02 1.88e-56 4.04e+53 1.44e+00 5.38e-03 8.21e-57 9.26e+53 + 2.00e+03 3 1.39e+01 2.10e-03 5.38e-03 7.70e+01 7.70e-05 1.10e-02 1.25e-02 3.34e+01 1.78e-04 2.70e-02 1.94e-02 7.46e+00 7.95e-04 4.99e-02 2.62e-02 8.73e-01 6.79e-03 7.94e-02 3.28e-02 5.46e-02 1.09e-01 1.15e-01 3.90e-02 1.86e-03 3.18e+00 1.57e-01 4.48e-02 3.59e-05 1.65e+02 2.05e-01 5.01e-02 4.03e-07 1.47e+04 2.57e-01 5.50e-02 2.74e-09 2.16e+06 3.14e-01 5.93e-02 1.18e-11 5.03e+08 3.76e-01 6.31e-02 3.37e-14 1.76e+11 4.40e-01 6.62e-02 6.75e-17 8.78e+13 5.08e-01 6.86e-02 1.00e-19 5.90e+16 5.77e-01 7.04e-02 1.19e-22 5.00e+19 6.48e-01 7.14e-02 1.19e-25 4.98e+22 7.20e-01 7.18e-02 1.09e-28 5.42e+25 7.92e-01 7.14e-02 9.88e-32 6.00e+28 8.63e-01 7.04e-02 9.47e-35 6.26e+31 9.32e-01 6.86e-02 1.03e-37 5.73e+34 1.00e+00 6.62e-02 1.38e-40 4.29e+37 1.06e+00 6.31e-02 2.42e-43 2.45e+40 1.13e+00 5.93e-02 5.91e-46 1.00e+43 1.18e+00 5.50e-02 2.14e-48 2.77e+45 1.24e+00 5.01e-02 1.22e-50 4.87e+47 1.28e+00 4.48e-02 1.14e-52 5.18e+49 1.32e+00 3.90e-02 1.86e-54 3.18e+51 1.36e+00 3.28e-02 5.48e-56 1.08e+53 1.39e+00 2.62e-02 3.01e-57 1.97e+54 1.41e+00 1.94e-02 3.18e-58 1.86e+55 1.43e+00 1.25e-02 6.62e-59 8.95e+55 1.44e+00 5.38e-03 2.75e-59 2.15e+56 + 3.00e+03 4 1.48e+01 2.10e-03 5.38e-03 7.71e+01 5.84e-05 1.10e-02 1.25e-02 3.34e+01 1.35e-04 2.70e-02 1.94e-02 7.44e+00 6.05e-04 4.99e-02 2.62e-02 8.68e-01 5.18e-03 7.94e-02 3.28e-02 5.39e-02 8.35e-02 1.15e-01 3.90e-02 1.82e-03 2.47e+00 1.57e-01 4.48e-02 3.45e-05 1.30e+02 2.05e-01 5.01e-02 3.76e-07 1.20e+04 2.57e-01 5.50e-02 2.44e-09 1.84e+06 3.14e-01 5.93e-02 9.75e-12 4.62e+08 3.76e-01 6.31e-02 2.48e-14 1.82e+11 4.40e-01 6.62e-02 4.19e-17 1.08e+14 5.08e-01 6.86e-02 4.95e-20 9.09e+16 5.77e-01 7.04e-02 4.39e-23 1.03e+20 6.48e-01 7.14e-02 3.16e-26 1.43e+23 7.20e-01 7.18e-02 2.01e-29 2.24e+26 7.92e-01 7.14e-02 1.24e-32 3.64e+29 8.63e-01 7.04e-02 8.00e-36 5.62e+32 9.32e-01 6.86e-02 5.91e-39 7.62e+35 1.00e+00 6.62e-02 5.38e-42 8.37e+38 1.06e+00 6.31e-02 6.49e-45 6.93e+41 1.13e+00 5.93e-02 1.11e-47 4.04e+44 1.18e+00 5.50e-02 2.90e-50 1.55e+47 1.24e+00 5.01e-02 1.22e-52 3.69e+49 1.28e+00 4.48e-02 8.71e-55 5.17e+51 1.32e+00 3.90e-02 1.11e-56 4.04e+53 1.36e+00 3.28e-02 2.66e-58 1.69e+55 1.39e+00 2.62e-02 1.23e-59 3.66e+56 1.41e+00 1.94e-02 1.14e-60 3.95e+57 1.43e+00 1.25e-02 2.16e-61 2.08e+58 1.44e+00 5.38e-03 8.54e-62 5.27e+58 + 4.00e+03 5 1.49e+01 2.10e-03 5.38e-03 7.86e+01 5.72e-05 1.10e-02 1.25e-02 3.34e+01 1.35e-04 2.70e-02 1.94e-02 7.16e+00 6.28e-04 4.99e-02 2.62e-02 7.83e-01 5.75e-03 7.94e-02 3.28e-02 4.39e-02 1.03e-01 1.15e-01 3.90e-02 1.27e-03 3.53e+00 1.57e-01 4.48e-02 1.94e-05 2.32e+02 2.05e-01 5.01e-02 1.58e-07 2.84e+04 2.57e-01 5.50e-02 7.21e-10 6.25e+06 3.14e-01 5.93e-02 1.93e-12 2.33e+09 3.76e-01 6.31e-02 3.28e-15 1.37e+12 4.40e-01 6.62e-02 3.78e-18 1.19e+15 5.08e-01 6.86e-02 3.21e-21 1.40e+18 5.77e-01 7.04e-02 2.16e-24 2.08e+21 6.48e-01 7.14e-02 1.25e-27 3.60e+24 7.20e-01 7.18e-02 6.67e-31 6.74e+27 7.92e-01 7.14e-02 3.56e-34 1.26e+31 8.63e-01 7.04e-02 2.05e-37 2.20e+34 9.32e-01 6.86e-02 1.36e-40 3.30e+37 1.00e+00 6.62e-02 1.13e-43 3.97e+40 1.06e+00 6.31e-02 1.26e-46 3.56e+43 1.13e+00 5.93e-02 2.02e-49 2.23e+46 1.18e+00 5.50e-02 4.92e-52 9.15e+48 1.24e+00 5.01e-02 1.94e-54 2.31e+51 1.28e+00 4.48e-02 1.32e-56 3.42e+53 1.32e+00 3.90e-02 1.61e-58 2.80e+55 1.36e+00 3.28e-02 3.69e-60 1.22e+57 1.39e+00 2.62e-02 1.66e-61 2.72e+58 1.41e+00 1.94e-02 1.50e-62 3.01e+59 1.43e+00 1.25e-02 2.79e-63 1.61e+60 1.44e+00 5.38e-03 1.09e-63 4.13e+60 diff --git a/src/fisst/COPYRIGHT b/src/fisst/COPYRIGHT new file mode 100644 index 0000000000..0a96a2eed9 --- /dev/null +++ b/src/fisst/COPYRIGHT @@ -0,0 +1,14 @@ +Copyright (c) 2020 of Glen Hocky + +The FISST module is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The FISST module is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with plumed. If not, see . diff --git a/src/fisst/FISST.cpp b/src/fisst/FISST.cpp new file mode 100644 index 0000000000..e5ec0682b9 --- /dev/null +++ b/src/fisst/FISST.cpp @@ -0,0 +1,661 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Copyright (c) 2020 of Glen Hocky + +The FISST module is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The FISST module is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#include "bias/Bias.h" +#include "core/ActionRegister.h" +#include "core/Atoms.h" +#include "core/PlumedMain.h" +#include "tools/File.h" +#include "tools/Matrix.h" +#include "tools/Random.h" +#include "legendre_rule_fast.h" + +#include + + +using namespace PLMD; +using namespace bias; + +//namespace is lowercase to match +//module names being all lowercase + +namespace PLMD { +namespace fisst { + +//+PLUMEDOC FISSTMOD_BIAS FISST +/* +Compute and apply the optimal linear force on an observable to enhance sampling +of conformational distributions over a range of applied forces. This method is described in \cite Hartmann-FISST-2019 + +If the system's Hamiltonian is given by: +\f[ + H(\vec{p},\vec{q}) = \sum_{j} \frac{p_j^2}{2m_j} + U(\vec{q}), +\f] + +This bias modifies the Hamiltonian to be: +\f[ + H'(\vec{p},\vec{q}) = H(\vec{p},\vec{q}) - \bar{F} Q +\f] + +where for CV \f$Q\f$, a coupling constant \f${\bar{F}}\f$ is determined +adaptively according to the FISST algorithm. + +Specifically, +\f[ +\bar{F}(Q)=\frac{ \int_{F_{min}}^{F_{max}} e^{\beta F Q(\vec{q})} \omega(F) F dF}{\int_{F_{min}}^{F_{max}} e^{\beta F Q(\vec{q})} \omega(F) dF}, +\f] + +where \f$\vec{q}\f$ are the molecular coordinates of the system, and \f$w(F)\f$ is a weighting function that is learned on the fly for each force by the FISST algorithm (starting from an initial weight distribution, uniform by default). + +The target for \f$w(F)=1/Z_q(F)\f$, where +\f[ + Z_q(F) \equiv \int d\vec{q} e^{-\beta U(\vec{q}) + \beta F Q(\vec{q})}. +\f] + +FISST also computes and writes Observable Weights \f$W_F(\vec{q}_t)\f$ for a molecular configuration at time \f$t\f$, so that averages of other quantities \f$A(\vec{q})\f$ can be reconstructed later at different force values (over a trajectory with \f$T\f$ samples): +\f[ + \langle A \rangle_F = \frac{1}{T} \sum_t W_F(\vec{q}_t) A(\vec{q}_t). +\f] + + +\par Examples + +In the following example, an adaptive restraint is learned to bias the distance between two atoms in a system, for a force range of 0-100 pN. + +\plumedfile +UNITS LENGTH=A TIME=fs ENERGY=kcal/mol + +b1: GROUP ATOMS=1 +b2: GROUP ATOMS=12 + +dend: DISTANCE ATOMS=b1,b2 + +#The conversion factor is 69.4786 pN = 1 kcal/mol/Angstrom + +#0 pN to 100 pN +f: FISST MIN_FORCE=0 MAX_FORCE=1.44 PERIOD=100 NINTERPOLATE=31 ARG=dend OUT_RESTART=pull.restart.txt OUT_OBSERVABLE=pull.observable.txt OBSERVABLE_FREQ=1000 + +PRINT ARG=dend,f.dend_fbar,f.bias,f.force2 FILE=pull.colvar.txt STRIDE=1000 +\endplumedfile + + +*/ +//+ENDPLUMEDOC + + +class FISST : public Bias { + + +private: + /*We will get this and store it once, since on-the-fly changing number of CVs will be fatal*/ + const unsigned int ncvs_; + std::vector center_; + std::vector current_avg_force_; + + std::vector forces_; + std::vector force_weight_; + std::vector gauss_weight_; + std::vector partition_estimate_; + std::vector observable_weight_; + + std::string in_restart_name_; + std::string out_restart_name_; + std::string out_observable_name_; + std::string fmt_; + std::string initial_weight_dist_; + OFile out_restart_; + OFile out_observable_; + IFile in_restart_; + bool b_freeze_; + bool b_adaptive_; + bool b_restart_; + bool b_write_restart_; + bool b_write_observable_; + bool b_first_restart_sample_; + int period_; + int reset_period_; + int observable_freq_; + int n_interpolation_; + int n_samples_; + double kbt_; + double beta_; + //change min_force and max_force to vectors if going to do more than one cv + double max_force_; + double min_force_; + double initial_weight_rate_; + double threshold_; + Random rand_; + + + Value* value_force2_; + void readInRestart(); + void NormalizeForceWeights(); + /*setup output restart*/ + void setupOutRestart(); + void setupOutObservable(); + /*write output restart*/ + void writeOutRestart(); + void writeOutObservable(); + void update_statistics(); + void update_bias(); + void apply_bias(); + void compute_observable_weight(); + +public: + explicit FISST(const ActionOptions&); + void calculate(); + void update(); + void turnOnDerivatives(); + static void registerKeywords(Keywords& keys); + ~FISST(); +}; + +PLUMED_REGISTER_ACTION(FISST,"FISST") + +void FISST::registerKeywords(Keywords& keys) { + Bias::registerKeywords(keys); + keys.use("ARG"); + keys.add("compulsory","PERIOD","Steps corresponding to the learning rate"); + keys.add("optional","RESET_PERIOD","Reset the learning statistics every time this number of steps comes around."); + keys.add("compulsory","NINTERPOLATE","Number of grid points on which to do interpolation."); + keys.add("compulsory","MIN_FORCE","Minimum force (per CV) to use for sampling. Units: [Energy]/[CV] (can be negative)."); + keys.add("compulsory","MAX_FORCE","Maximum force (per CV) to use for sampling."); + keys.add("compulsory","CENTER","0","The CV value at which the applied bias energy will be zero"); + keys.add("optional","KBT","The system temperature in units of KB*T. If not provided will be taken from MD code (if available)"); + + keys.add("optional","INITIAL_WEIGHT_DIST","Starting distribution for the force weights (options: UNIFORM, EXP, GAUSS)."); + keys.add("optional","INITIAL_WEIGHT_RATE","Rate of decay for exponential and gaussian distributions. W(F)~exp(-r |F|^d)."); + + keys.add("optional","RESTART_FMT","the format that should be used to output real numbers in FISST restarts."); + keys.add("optional","OUT_RESTART","Output file for all information needed to continue FISST simulation." + "If you have the RESTART directive set (global or for FISST), this file will be appended to." + "Note that the header will be printed again if appending."); + keys.add("optional","IN_RESTART","Read this file to continue an FISST simulation. " + "If same as OUT_RESTART and you have not set the RESTART directive, the file will be backed-up and overwritten with new output." + "If you do have the RESTART flag set and it is the same name as OUT_RESTART, this file will be appended."); + keys.add("optional","OUT_OBSERVABLE","Output file putting weights needed to compute observables at different force values." + "If you have the RESTART directive set (global or for FISST), this file will be appended to. " + "Note that the header will be printed again if appending."); + keys.add("optional","OBSERVABLE_FREQ","How often to write out observable weights (default=period)."); + keys.addFlag("FREEZE",false,"Fix bias weights at current level (only used for restarting)."); + keys.use("RESTART"); + keys.addOutputComponent("force2","default","squared value of force from the bias."); + keys.addOutputComponent("_fbar","default", "For each named CV biased, there will be a corresponding output CV_fbar storing the current linear bias prefactor."); +} + +FISST::FISST(const ActionOptions&ao): + PLUMED_BIAS_INIT(ao), + ncvs_(getNumberOfArguments()), + current_avg_force_(ncvs_,0.0), + center_(ncvs_,0.0), + //change min_force and max_force to vectors if going to do more than one cv + min_force_(0.0), + max_force_(0.0), + in_restart_name_(""), + out_restart_name_(""), + out_observable_name_(""), + fmt_("%e"), + b_freeze_(false), + b_restart_(false), + b_write_restart_(false), + b_write_observable_(false), + b_first_restart_sample_(true), + n_interpolation_(0), + n_samples_(0), + initial_weight_rate_(0), + initial_weight_dist_("UNIFORM"), + period_(0), + reset_period_(0), + observable_freq_(0), + kbt_(0.0), + value_force2_(NULL) +{ + if(ncvs_==0) + error("Must specify at least one CV with ARG"); + + //temporary + if(ncvs_>1) + error("FISST only supports using one CV right now"); + + addComponent("force2"); + componentIsNotPeriodic("force2"); + value_force2_ = getPntrToComponent("force2"); + + for(unsigned int i = 0; igetName() + "_fbar"; + addComponent(comp); + componentIsNotPeriodic(comp); + } + + parseVector("CENTER",center_); + //change min_force and max_force to vectors if going to do more than one cv + parse("MIN_FORCE",min_force_); + parse("MAX_FORCE",max_force_); + parse("PERIOD",period_); + parse("RESET_PERIOD",reset_period_); + parse("INITIAL_WEIGHT_DIST",initial_weight_dist_); + parse("INITIAL_WEIGHT_RATE",initial_weight_rate_); + parse("OBSERVABLE_FREQ",observable_freq_); + parse("NINTERPOLATE",n_interpolation_); + parseFlag("FREEZE",b_freeze_); + parse("KBT",kbt_); + parse("RESTART_FMT", fmt_); + fmt_ = " " + fmt_;//add space since parse strips them + parse("OUT_RESTART",out_restart_name_); + parse("OUT_OBSERVABLE",out_observable_name_); + parse("IN_RESTART",in_restart_name_); + checkRead(); + + if(center_.size() != ncvs_) + error("Must have same number of CENTER arguments as ARG arguments"); + + if(in_restart_name_ != "") { + b_restart_ = true; + log.printf(" reading simulation information from file: %s\n",in_restart_name_.c_str()); + readInRestart(); + } else { + + if(! kbt_ > 0.0) + kbt_ = plumed.getAtoms().getKbT(); + + //in driver, this results in kbt of 0 + if(kbt_ == 0) { + error(" Unable to determine valid kBT. " + "Could be because you are runnning from driver or MD didn't give temperature.\n" + "Consider setting temperature manually with the KBT keyword."); + } + + log.printf(" kBT = %f\n",kbt_); + log.printf(" Updating with a time scale of %i steps\n",period_); + + log.printf(" Using centers for CVs of:"); + for(unsigned int i = 0; i< ncvs_; i++) { + log.printf(" %f ",center_[i]); + } + log.printf("\n"); + observable_weight_.resize(n_interpolation_); + for(unsigned int i = 0; i0) { + log.printf(" writing restart information every %i steps to file %s with format %s\n",abs(period_),out_restart_name_.c_str(), fmt_.c_str()); + b_write_restart_ = true; + setupOutRestart(); + } + if(out_observable_name_.length()>0) { + if(observable_freq_==0) observable_freq_ = period_; + log.printf(" writing observable information every %i steps to file %s with format %s\n",observable_freq_,out_observable_name_.c_str(), fmt_.c_str()); + b_write_observable_ = true; + setupOutObservable(); + } + + //add citation later: + //log<<" Bibliography "<getName(); + in_restart_.scanField(cv_name,tmp); + for(unsigned int j =0; j0) out_restart_.addConstantField("reset_period").printField("reset_period",reset_period_); + out_restart_.addConstantField("min_force").printField("min_force",min_force_); + out_restart_.addConstantField("max_force").printField("max_force",max_force_); +} + +void FISST::writeOutRestart() { + std::string cv_name; + out_restart_.printField("time",getTimeStep()*getStep()); + out_restart_.printField("nsamples",n_samples_); + + for(unsigned int i = 0; igetName(); + double Q_i = difference(i, center_[i], getArgument(i)); + out_restart_.printField(cv_name,Q_i); + for(int j = 0; j < n_interpolation_; j++ ) { +//have to update this for multiple cvs + out_restart_.printField(cv_name + "_f"+std::to_string(j),forces_[j]); + out_restart_.printField(cv_name + "_g"+std::to_string(j),gauss_weight_[j]); + out_restart_.printField(cv_name + "_w"+std::to_string(j),force_weight_[j]); + out_restart_.printField(cv_name + "_z"+std::to_string(j),partition_estimate_[j]); + } + } + out_restart_.printField(); +} + +void FISST::writeOutObservable() { + std::string cv_name; + out_observable_.printField("time",getTimeStep()*getStep()); + out_observable_.printField("nsamples",n_samples_); + + for(unsigned int i = 0; igetName(); + double Q_i = difference(i, center_[i], getArgument(i)); + out_observable_.printField(cv_name,Q_i); + out_observable_.printField(cv_name + "_fbar",current_avg_force_[i]); + for(int j = 0; j < n_interpolation_; j++ ) { +//have to update this for multiple cvs + out_observable_.printField(cv_name + "_f"+std::to_string(j),forces_[j]); + out_observable_.printField(cv_name + "_ow"+std::to_string(j),observable_weight_[j]); + } + } + out_observable_.printField(); +} + + +void FISST::calculate() { + if(getStep() == 0 ) { + if(b_write_restart_) writeOutRestart(); + if(b_write_observable_) writeOutObservable(); + } + + if(! b_freeze_) { + if(b_restart_ && b_first_restart_sample_) { + //dont' update statistics if restarting and first sample + b_first_restart_sample_ = false; + } + else { + update_statistics(); + } + } + update_bias(); + apply_bias(); + + //check about writing restart file + if(getStep()>0 && getStep()%period_==0) { + if(b_write_restart_) writeOutRestart(); + } + if(getStep()>0 && getStep()%observable_freq_==0) { + if(b_write_observable_) { + compute_observable_weight(); + writeOutObservable(); + } + } +} + + +void FISST::apply_bias() { + //Compute linear force as in "restraint" + double ene = 0, totf2 = 0, cv, m, f; + + for(unsigned int i = 0; i < ncvs_; ++i) { + cv = difference(i, center_[i], getArgument(i)); + double fbar = current_avg_force_[i]; + ene -= fbar*cv; + setOutputForce(i,fbar); + totf2 += fbar*fbar; + + std::string fbar_name_ = getPntrToArgument(i)->getName() + "_fbar"; + Value* fbar_ = getPntrToComponent(fbar_name_); + fbar_->set(fbar); + }; + + setBias(ene); + value_force2_->set(totf2); + //log.flush(); +} + +void FISST::update_statistics() { +//get stride is for multiple time stepping + double dt=getTimeStep()*getStride(); + double h = dt/(period_*getTimeStep()); + double fbar_denum_integral = 0.0; + + int step = getStep(); + if(reset_period_>0 && step>0 && step%reset_period_==0) { + n_samples_=1; + } + else { + n_samples_++; + } + double d_n_samples = (double)n_samples_; + + for(unsigned int i = 0; i < ncvs_; ++i) { + double Q_i = difference(i, center_[i], getArgument(i)); + for(unsigned int j=0; j. ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +//#ifndef __PLUMED_fisst_legendre_rule_fast_h +//#define __PLUMED_fisst_legendre_rule_fast_h +// adapted from: +// https://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule_fast/legendre_rule_fast.html +// +#include "legendre_rule_fast.h" + +namespace PLMD { +namespace fisst { +//****************************************************************************80 + +void legendre_compute_glr ( int n, double x[], double w[] ) + +//****************************************************************************80 +// +// Purpose: +// +// LEGENDRE_COMPUTE_GLR: Legendre quadrature by the Glaser-Liu-Rokhlin method. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 20 October 2009 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Reference: +// +// Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin, +// A fast algorithm for the calculation of the roots of special functions, +// SIAM Journal on Scientific Computing, +// Volume 29, Number 4, pages 1420-1438, 2007. +// +// Parameters: +// +// Input, int N, the order. +// +// Output, double X[N], the abscissas. +// +// Output, double W[N], the weights. +// +{ + int i; + double p; + double pp; + double w_sum; +// +// Get the value and derivative of the N-th Legendre polynomial at 0. +// + legendre_compute_glr0 ( n, &p, &pp ); +// +// If N is odd, then zero is a root. +// + if ( n % 2 == 1 ) + { + x[(n-1)/2] = p; + w[(n-1)/2] = pp; + } +// +// If N is even, we have to call a function to find the first root. +// + else + { + legendre_compute_glr2 ( p, n, &x[n/2], &w[n/2] ); + } +// +// Get the complete set of roots and derivatives. +// + legendre_compute_glr1 ( n, x, w ); +// +// Compute the W. +// + for ( i = 0; i < n; i++ ) + { + w[i] = 2.0 / ( 1.0 - x[i] ) / ( 1.0 + x[i] ) / w[i] / w[i]; + } + w_sum = 0.0; + for ( i = 0; i < n; i++ ) + { + w_sum = w_sum + w[i]; + } + for ( i = 0; i < n; i++ ) + { + w[i] = 2.0 * w[i] / w_sum; + } + return; +} +//****************************************************************************80 + +void legendre_compute_glr0 ( int n, double *p, double *pp ) + +//****************************************************************************80 +// +// Purpose: +// +// LEGENDRE_COMPUTE_GLR0 gets a starting value for the fast algorithm. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 19 October 2009 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Reference: +// +// Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin, +// A fast algorithm for the calculation of the roots of special functions, +// SIAM Journal on Scientific Computing, +// Volume 29, Number 4, pages 1420-1438, 2007. +// +// Parameters: +// +// Input, int N, the order of the Legendre polynomial. +// +// Output, double *P, *PP, the value of the N-th Legendre polynomial +// and its derivative at 0. +// +{ + double dk; + int k; + double pm1; + double pm2; + double ppm1; + double ppm2; + + pm2 = 0.0; + pm1 = 1.0; + ppm2 = 0.0; + ppm1 = 0.0; + + for ( k = 0; k < n; k++) + { + dk = ( double ) k; + *p = - dk * pm2 / ( dk + 1.0 ); + *pp = ( ( 2.0 * dk + 1.0 ) * pm1 - dk * ppm2 ) / ( dk + 1.0 ); + pm2 = pm1; + pm1 = *p; + ppm2 = ppm1; + ppm1 = *pp; + } + return; +} +//****************************************************************************80 + +void legendre_compute_glr1 ( int n, double *x, double *w ) + +//****************************************************************************80 +// +// Purpose: +// +// LEGENDRE_COMPUTE_GLR1 gets the complete set of Legendre points and weights. +// +// Discussion: +// +// This routine requires that a starting estimate be provided for one +// root and its derivative. This information will be stored in entry +// (N+1)/2 if N is odd, or N/2 if N is even, of X and W. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 19 October 2009 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Reference: +// +// Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin, +// A fast algorithm for the calculation of the roots of special functions, +// SIAM Journal on Scientific Computing, +// Volume 29, Number 4, pages 1420-1438, 2007. +// +// Parameters: +// +// Input, int N, the order of the Legendre polynomial. +// +// Input/output, double X[N]. On input, a starting value +// has been set in one entry. On output, the roots of the Legendre +// polynomial. +// +// Input/output, double W[N]. On input, a starting value +// has been set in one entry. On output, the derivatives of the Legendre +// polynomial at the zeros. +// +// Local Parameters: +// +// Local, int M, the number of terms in the Taylor expansion. +// +{ + double dk; + double dn; + double h; + int j; + int k; + int l; + int m = 30; + int n2; + static double pi = 3.141592653589793; + int s; + double *u; + double *up; + double xp; + + if ( n % 2 == 1 ) + { + n2 = ( n - 1 ) / 2 - 1; + s = 1; + } + else + { + n2 = n / 2 - 1; + s = 0; + } + + u = new double[m+2]; + up = new double[m+1]; + + dn = ( double ) n; + + for ( j = n2 + 1; j < n - 1; j++ ) + { + xp = x[j]; + + h = rk2_leg ( pi/2.0, -pi/2.0, xp, n ) - xp; + + u[0] = 0.0; + u[1] = 0.0; + u[2] = w[j]; + + up[0] = 0.0; + up[1] = u[2]; + + for ( k = 0; k <= m - 2; k++ ) + { + dk = ( double ) k; + + u[k+3] = + ( + 2.0 * xp * ( dk + 1.0 ) * u[k+2] + + ( dk * ( dk + 1.0 ) - dn * ( dn + 1.0 ) ) * u[k+1] / ( dk + 1.0 ) + ) / ( 1.0 - xp ) / ( 1.0 + xp ) / ( dk + 2.0 ); + + up[k+2] = ( dk + 2.0 ) * u[k+3]; + } + + for ( l = 0; l < 5; l++ ) + { + h = h - ts_mult ( u, h, m ) / ts_mult ( up, h, m-1 ); + } + + x[j+1] = xp + h; + w[j+1] = ts_mult ( up, h, m - 1 ); + } + + for ( k = 0; k <= n2 + s; k++ ) + { + x[k] = - x[n-1-k]; + w[k] = w[n-1-k]; + } + return; +} +//****************************************************************************80 + +void legendre_compute_glr2 ( double pn0, int n, double *x1, double *d1 ) + +//****************************************************************************80 +// +// Purpose: +// +// LEGENDRE_COMPUTE_GLR2 finds the first real root. +// +// Discussion: +// +// This function is only called if N is even. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 19 October 2009 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Reference: +// +// Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin, +// A fast algorithm for the calculation of the roots of special functions, +// SIAM Journal on Scientific Computing, +// Volume 29, Number 4, pages 1420-1438, 2007. +// +// Parameters: +// +// Input, double PN0, the value of the N-th Legendre polynomial +// at 0. +// +// Input, int N, the order of the Legendre polynomial. +// +// Output, double *X1, the first real root. +// +// Output, double *D1, the derivative at X1. +// +// Local Parameters: +// +// Local, int M, the number of terms in the Taylor expansion. +// +{ + double dk; + double dn; + int k; + int l; + int m = 30; + static double pi = 3.141592653589793; + double t; + double *u; + double *up; + + t = 0.0; + *x1 = rk2_leg ( t, -pi/2.0, 0.0, n ); + + u = new double[m+2]; + up = new double[m+1]; + + dn = ( double ) n; +// +// U[0] and UP[0] are never used. +// U[M+1] is set, but not used, and UP[M] is set and not used. +// What gives? +// + u[0] = 0.0; + u[1] = pn0; + + up[0] = 0.0; + + for ( k = 0; k <= m - 2; k = k + 2 ) + { + dk = ( double ) k; + + u[k+2] = 0.0; + u[k+3] = ( dk * ( dk + 1.0 ) - dn * ( dn + 1.0 ) ) * u[k+1] + / (dk + 1.0) / (dk + 2.0 ); + + up[k+1] = 0.0; + up[k+2] = ( dk + 2.0 ) * u[k+3]; + } + + for ( l = 0; l < 5; l++ ) + { + *x1 = *x1 - ts_mult ( u, *x1, m ) / ts_mult ( up, *x1, m-1 ); + } + *d1 = ts_mult ( up, *x1, m-1 ); + + return; +} + +//****************************************************************************80 + +void rescale ( double a, double b, int n, double x[], double w[] ) + +//****************************************************************************80 +// +// Purpose: +// +// RESCALE rescales a Legendre quadrature rule from [-1,+1] to [A,B]. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 18 October 2009 +// +// Author: +// +// Original MATLAB version by Nick Hale. +// C++ version by John Burkardt. +// +// Reference: +// +// Andreas Glaser, Xiangtao Liu, Vladimir Rokhlin, +// A fast algorithm for the calculation of the roots of special functions, +// SIAM Journal on Scientific Computing, +// Volume 29, Number 4, pages 1420-1438, 2007. +// +// Parameters: +// +// Input, double A, B, the endpoints of the new interval. +// +// Input, int N, the order. +// +// Input/output, double X[N], on input, the abscissas for [-1,+1]. +// On output, the abscissas for [A,B]. +// +// Input/output, double W[N], on input, the weights for [-1,+1]. +// On output, the weights for [A,B]. +// +{ + int i; + + for ( i = 0; i < n; i++ ) + { + x[i] = ( ( a + b ) + ( b - a ) * x[i] ) / 2.0; + } + for ( i = 0; i < n; i++ ) + { + w[i] = ( b - a ) * w[i] / 2.0; + } + return; +} +//****************************************************************************80 + +double rk2_leg ( double t1, double t2, double x, int n ) + +//****************************************************************************80 +// +// Purpose: +// +// RK2_LEG advances the value of X(T) using a Runge-Kutta method. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 22 October 2009 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Parameters: +// +// Input, double T1, T2, the range of the integration interval. +// +// Input, double X, the value of X at T1. +// +// Input, int N, the number of steps to take. +// +// Output, double RK2_LEG, the value of X at T2. +// +{ + double f; + double h; + int j; + double k1; + double k2; + int m = 10; + double snn1; + double t; + + h = ( t2 - t1 ) / ( double ) m; + snn1 = sqrt ( ( double ) ( n * ( n + 1 ) ) ); + t = t1; + + for ( j = 0; j < m; j++ ) + { + f = ( 1.0 - x ) * ( 1.0 + x ); + k1 = - h * f / ( snn1 * sqrt ( f ) - 0.5 * x * sin ( 2.0 * t ) ); + x = x + k1; + + t = t + h; + + f = ( 1.0 - x ) * ( 1.0 + x ); + k2 = - h * f / ( snn1 * sqrt ( f ) - 0.5 * x * sin ( 2.0 * t ) ); + x = x + 0.5 * ( k2 - k1 ); + } + return x; +} +//****************************************************************************80 + +double ts_mult ( double *u, double h, int n ) + +//****************************************************************************80 +// +// Purpose: +// +// TS_MULT evaluates a polynomial. +// +// Licensing: +// +// This code is distributed under the GNU LGPL license. +// +// Modified: +// +// 17 May 2013 +// +// Author: +// +// Original C++ version by Nick Hale. +// This C++ version by John Burkardt. +// +// Parameters: +// +// Input, double U[N+1], the polynomial coefficients. +// U[0] is ignored. +// +// Input, double H, the polynomial argument. +// +// Input, int N, the number of terms to compute. +// +// Output, double TS_MULT, the value of the polynomial. +// +{ + double hk; + int k; + double ts; + + ts = 0.0; + hk = 1.0; + for ( k = 1; k<= n; k++ ) + { + ts = ts + u[k] * hk; + hk = hk * h; + } + return ts; +} + +//close the namespaces +} +} + +//#endif + diff --git a/src/fisst/legendre_rule_fast.h b/src/fisst/legendre_rule_fast.h new file mode 100644 index 0000000000..99433c010d --- /dev/null +++ b/src/fisst/legendre_rule_fast.h @@ -0,0 +1,44 @@ +/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +Copyright (c) 2020 of Glen Hocky + +The FISST module is free software: you can redistribute it and/or modify +it under the terms of the GNU Lesser General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The FISST module is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU Lesser General Public License for more details. + +You should have received a copy of the GNU Lesser General Public License +along with plumed. If not, see . ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ +#ifndef __PLUMED_fisst_legendre_rule_fast_h +#define __PLUMED_fisst_legendre_rule_fast_h +// taken from: +// https://people.sc.fsu.edu/~jburkardt/cpp_src/legendre_rule_fast/legendre_rule_fast.html +// +# include +# include +# include +# include +# include +# include +# include + +namespace PLMD { +namespace fisst { + +void legendre_compute_glr ( int n, double x[], double w[] ); +void legendre_compute_glr0 ( int n, double *p, double *pp ); +void legendre_compute_glr1 ( int n, double *roots, double *ders ); +void legendre_compute_glr2 ( double p, int n, double *roots, double *ders ); +void rescale ( double a, double b, int n, double x[], double w[] ); +double rk2_leg ( double t, double tn, double x, int n ); +double ts_mult ( double *u, double h, int n ); + +} +} + +#endif diff --git a/src/fisst/module.type b/src/fisst/module.type new file mode 100644 index 0000000000..de83273033 --- /dev/null +++ b/src/fisst/module.type @@ -0,0 +1 @@ +default-off diff --git a/user-doc/FISSTMOD.md b/user-doc/FISSTMOD.md new file mode 100644 index 0000000000..0921e23a61 --- /dev/null +++ b/user-doc/FISSTMOD.md @@ -0,0 +1,26 @@ +\page FISSTMOD FISST (Infinite Switch Simulated Tempering in Force) + + + +## Overview + +This FISST module contains methods for adaptively determining weight parameters to construct a bias function that represents the Infinite Switch limit of Simulated Tempering for a linear bias coefficient of a CV, as described in \cite Hartmann-FISST-2019. + +## Installation +This module is not installed by default. Add '\-\-enable-modules=fisst' to your './configure' command when building PLUMED to enable these features. + +## Usage +Currently, all features of the FISST module are included in a single FISST bias function: \ref FISST + +## Module Contents +- \subpage FISSTMODBias + +\page FISSTMODBias Biases Documentation + +The following list contains descriptions of biases developed for the PLUMED-FISST module. They can be used in combination with other biases outside of the FISST module. + +@FISSTMOD_BIAS@ diff --git a/user-doc/bibliography.bib b/user-doc/bibliography.bib index 47b36bbc07..c7d50749bd 100644 --- a/user-doc/bibliography.bib +++ b/user-doc/bibliography.bib @@ -2990,38 +2990,43 @@ @article{bussi2015free } @article{Giberti_jp16100, -author = {Giberti, F. and Cheng, B. and Tribello, G. A. and Ceriotti, M.}, -title = {Iterative Unbiasing of Quasi-Equilibrium Sampling}, -journal = {Journal of Chemical Theory and Computation}, -volume = {16}, -number = {1}, -pages = {100-107}, -year = {2020}, -doi = {10.1021/acs.jctc.9b00907}, - note ={PMID: 31743021}, + author = {Giberti, F. and Cheng, B. and Tribello, G. A. and Ceriotti, M.}, + title = {Iterative Unbiasing of Quasi-Equilibrium Sampling}, + journal = {Journal of Chemical Theory and Computation}, + volume = {16}, + number = {1}, + pages = {100-107}, + year = {2020}, + doi = {10.1021/acs.jctc.9b00907}, + note ={PMID: 31743021}, + URL = { + https://doi.org/10.1021/acs.jctc.9b00907 -URL = { - https://doi.org/10.1021/acs.jctc.9b00907 - -}, -eprint = { - https://doi.org/10.1021/acs.jctc.9b00907 - -} + }, + eprint = { + https://doi.org/10.1021/acs.jctc.9b00907 + } } @article{Hovan2019, -author = {Hovan, Ladislav and Comitani, Federico and Gervasio, Francesco L}, -doi = {10.1021/acs.jctc.8b00563}, -issn = {1549-9618}, -journal = {Journal of Chemical Theory and Computation}, -number = {1}, -pages = {25--32}, -publisher = {American Chemical Society}, -title = {{An Optimal Metric for the Path Collective Variables}}, -volume = {15}, -year = {2019} + author = {Hovan, Ladislav and Comitani, Federico and Gervasio, Francesco L}, + doi = {10.1021/acs.jctc.8b00563}, + issn = {1549-9618}, + journal = {Journal of Chemical Theory and Computation}, + number = {1}, + pages = {25--32}, + publisher = {American Chemical Society}, + title = {{An Optimal Metric for the Path Collective Variables}}, + volume = {15}, + year = {2019} +} + +@article{Hartmann-FISST-2019, + title={Infinite Switch Simulated Tempering in Force (FISST)}, + author={Hartmann, Michael J and Singh, Yuvraj and Vanden-Eijnden, Eric and Hocky, Glen M}, + journal={arXiv:1910.14064}, + year={2019}, } @Comment{jabref-meta: databaseType:bibtex;}