From 503cab3258bbe63881a655974692708b93dc240e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=A9r=C3=B4me=20H=C3=A9nin?= Date: Thu, 18 Jan 2024 15:15:49 +0100 Subject: [PATCH] Squashed commit of the following: MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit commit 0573ac25a8478c00704302d7ab5bb132c2352c99 Merge: 2b3f2a2c c2bdeb95 Author: Jérôme Hénin Date: Wed Jan 17 16:11:14 2024 +0100 Merge branch 'master' into new_abf_for_merging commit 2b3f2a2cb4d9c67e5c686428e4a7c7a9d0a62457 Author: Jérôme Hénin Date: Wed Jan 17 16:03:57 2024 +0100 Remove exptl feature linABF for now Needs to be better tested and documented commit fc86e9a94072b32ace0d1f8d5a0aa459859e7416 Author: Jérôme Hénin Date: Fri Jan 5 17:57:00 2024 +0100 Fix memory issues found by ASAN - alloc-dealloc mismatch - missing deallocation of features in derived class colvarbias_abf commit f88563846c989df5536b8be3129a3f98bcddf331 Author: Jérôme Hénin Date: Fri Jan 5 17:42:45 2024 +0100 Amend mwABF test script for NAMD-netlrts commit e1bf0c0dc64d247d1115a3bab13d5c4c0790f696 Author: Jérôme Hénin Date: Fri Jan 5 17:12:07 2024 +0100 Update mwABF documentation commit c88525b9445eac1aed6334a23ee781d2a3949131 Merge: 15a40973 6bb881bb Author: Jérôme Hénin Date: Wed Nov 29 12:15:33 2023 +0100 Merge branch 'master' into new_abf_for_merging commit 15a4097352f442430f57edc92799bd99c70e5ec5 Author: Jérôme Hénin Date: Wed Nov 1 21:39:39 2023 +0100 Adjust regtest for correct behavior of new ABF code as in d4b2b4cc1d7b35d87640824f3998eacd027038d0 commit 0aca6e8b2d410e1ef7ea26ababa7aaaa47c5a576 Author: Jérôme Hénin Date: Tue Oct 31 19:06:54 2023 +0100 Appease compiler warnings - "possible" memory leak (really not possible, but what can I do) - "garbage values" which weren't but it took some persuading - hidden variable (that was actually a nice catch by the Sun compiler) commit 0524fcc7b8ecbf2e3eb62dbb3c6330c8e73c9761 Author: Jérôme Hénin Date: Tue Oct 31 17:52:04 2023 +0100 Add missing override qualifiers commit 5a2c41a01ce06cb5da0633e3fc54ab19976760d1 Author: Jérôme Hénin Date: Mon Oct 30 15:27:37 2023 +0100 Remove branch README commit 2e37866330c6adbd27eb87cb7094701ac8c8d45b Author: Jérôme Hénin Date: Mon Oct 30 15:24:19 2023 +0100 Add include for smart pointers commit 0e3ff84a8687d2b674e0d4d775537dd4f0183aff Author: Jérôme Hénin Date: Mon Oct 30 11:38:31 2023 +0100 Remove redundant lines from merge commit 26cf5e8d85470ffe282595d4513945d8cbc1b69a Author: Jérôme Hénin Date: Thu Oct 26 17:27:55 2023 +0200 Remove dead code from ABF: adaptation_time, weights commit d4b2b4cc1d7b35d87640824f3998eacd027038d0 Author: Jérôme Hénin Date: Thu Oct 26 16:46:59 2023 +0200 Adjust regtest reference for new ABF code The new code has the correct behavior, which is to ignore the duplicate total force value obtained at the beginning of a second run statement. commit 6cb1f19e6bb2b708698b402fcf589215ab9cc964 Author: Jérôme Hénin Date: Thu Oct 26 14:50:49 2023 +0200 Remove unnecessary include commit ee34329774dd6fd224c7f1c7bf98eea545c5887a Author: Jérôme Hénin Date: Wed Oct 18 17:19:49 2023 +0200 Fix typos commit fdc5fb551f908f37b32e54c0b92f5215ae9c290e Author: Jérôme Hénin Date: Wed Oct 18 17:15:09 2023 +0200 Merge changes to NAMD tests from master commit 268f63c8ad48586f4c08d7dad6ed46715dae9e66 Author: Jérôme Hénin Date: Fri Oct 13 16:58:57 2023 +0200 Remove experimental ABF-AR features commit 32719e52f3e33438843c576ef8a6f1339f67b020 Merge: 1d5c2d97 46f1d6d7 Author: Jérôme Hénin Date: Fri Oct 13 11:30:16 2023 +0200 Merge branch 'master' into new_abf_for_merging commit 1d5c2d97333f90e80301b0808fd4d36195909b59 Merge: e396b025 7318e2ce Author: Jérôme Hénin Date: Fri Oct 6 19:03:12 2023 +0200 Merge branch 'master' into new_abf_for_merging commit e396b02566e6d04c37f47375172299343fa29274 Merge: 5e01a4f9 ab60dcb6 Author: Jérôme Hénin Date: Fri Oct 6 16:38:57 2023 +0200 Merge branch 'master' into new_abf_for_merging commit 5e01a4f9cf458ca42c200448160e8d478e34d501 Author: Jérôme Hénin Date: Fri Nov 11 23:55:02 2022 +0100 Make ABF grids unique_ptr commit 8efc8f6a0e73759b591e93fdac01d4f43a5521f1 Merge: 5e155e6e 3af3d4ed Author: Jérôme Hénin Date: Fri Nov 11 21:05:47 2022 +0100 Merge branch 'master' into new_abf commit 5e155e6e217e23102274e232fc746102dc591f0f Author: Jérôme Hénin Date: Fri Aug 12 23:28:59 2022 +0200 Initial version of README for branch changes commit 3bb83e2c38a3beb020146a35facfe54d3a6a6275 Author: Jérôme Hénin Date: Fri Aug 12 23:17:24 2022 +0200 Remove uninitialized b_linabf Boolean commit 19a746464d0bace9a090908d1abd0c030ad14bf7 Merge: b215af11 732a460d Author: Jérôme Hénin Date: Fri Aug 12 22:17:35 2022 +0200 Merge branch 'master' into new_abf commit b215af112297b6d9aa337693130e686b033b00cc Author: Jérôme Hénin Date: Tue Mar 29 18:14:00 2022 +0200 Add one more fragment lost in merge (alternate integrate potential constructor) commit 44dc35c564bedb80027c1e9efc5bd1a4372f309a Author: Jérôme Hénin Date: Tue Mar 29 18:05:00 2022 +0200 Update script doc commit d0cb0ba641161e41f94753afef62c129d6b1fc46 Merge: 7770165a 3d6c0e63 Author: Jérôme Hénin Date: Tue Mar 29 17:38:17 2022 +0200 Merge branch 'master' into new_abf commit 7770165aa0f950c3c8b30694d634fd5c7640f2c6 Author: Jérôme Hénin Date: Mon Mar 28 14:37:35 2022 +0200 Add fragment lost in merge commit 37b4124ec48d8d7f79101d445bed6402d598ea74 Author: Jérôme Hénin Date: Thu Mar 24 19:02:54 2022 +0100 Fix error introduced by merge commit 6255aa42f1a149b9ec8f68dfb6d10dcf70d0238f Author: Jérôme Hénin Date: Thu Mar 24 18:54:33 2022 +0100 Fix compiler warnings commit 420eda8f8811e9ba0708f323ea14fc9fb7184188 Merge: ac0d1896 f2a0dff6 Author: Jérôme Hénin Date: Thu Mar 24 18:31:52 2022 +0100 Merge branch 'master' into new_abf commit ac0d18966adcbf7ec21bd6bda3da291dd57066ef Author: Jérôme Hénin Date: Mon Jun 28 14:44:15 2021 +0200 WIP on ABF-AX: communicate and store local data commit f0c836fa9d88b3bb48a0fd9c93284bd27964ba7e Author: Jérôme Hénin Date: Wed Jun 23 14:44:10 2021 +0200 WIP on mwABF-AX commit 115f946e02412d3e26aa0419cb834b0a3d3beb99 Author: Jérôme Hénin Date: Sun Dec 20 20:57:51 2020 +0100 First version of Adiabatic Exploration commit e43c8aa125d08b92cb1f6c23cf8b27b46a2f658d Author: Jérôme Hénin Date: Fri Nov 6 11:40:58 2020 +0100 Parse updateBias after ABF-AR commit 2df66d3a69a5357377f116ce3b2a3772c54da5dd Author: Jérôme Hénin Date: Thu Nov 5 12:42:26 2020 +0100 WIP commit 6aa102c0498acec7fcc273c17dcebeb8992e6609 Author: Jérôme Hénin Date: Tue Oct 20 11:51:24 2020 +0200 ABF-AR in dimension > 1 Including lambda contributions beyond grid edges. also: - user control of minSamples - still playing with adaptation parameter of exploration bias - calculate ABF bias energy beyond grid edges commit 989978768dc92a555cd7cf7daead4beee69660f8 Author: Jérôme Hénin Date: Sat Oct 17 21:45:57 2020 +0200 First step towards multidim ABF-AR commit 4954706330377482b3e7ecddabffb475236ec689 Author: Jérôme Hénin Date: Fri Oct 16 13:12:39 2020 +0200 Reset ext coordinate after jumps in colvar this would in particular break mw-eABF with selection as the selection process creates jumps in Cartesian coords and hence in colvar values. commit cb1c0f102ed69f2a5e0ad86a6aa79a46a4088fb3 Author: Jérôme Hénin Date: Thu Oct 8 17:00:17 2020 +0200 Common smoothing functions for ABF and ABF-AR using either fullSamples or adaptationTime Also implemented choice of cosine or Gaussian kernel in ABF-AR commit 641242b4d34536436202550abc7f175f9b655240 Author: Jérôme Hénin Date: Thu Oct 1 12:27:53 2020 +0200 WIP on ABF-AR - handle special file names in colvars_grid.py - introduce cosine kernel - rename adaptationRate to adaptationTime commit 57f7f25eb99046ac236bc9999bda2ed23d688659 Author: Jérôme Hénin Date: Tue Sep 29 09:59:14 2020 +0200 WIP on optimization of ABF-AR loop, BC commit 79d7e6f00dc1be45ddbb25ccc4713e1cb483b7a9 Author: Jérôme Hénin Date: Sat Sep 26 17:46:44 2020 +0200 Update NAMD ABF test ref output for proper handling of zero step data previous reference had redundant data for one step commit c5ac01ce27a47fc6040d01ebe1c67942ee14b8b5 Author: Jérôme Hénin Date: Sat Sep 26 17:42:50 2020 +0200 Enable dependency tree items for ABF more generally, for derived classes of colvarbias commit f871fed5aca6094b4c9a75d204ba05fdaa69f59d Author: Jérôme Hénin Date: Fri Sep 25 18:53:19 2020 +0200 Add back missing changes after rebase commit 9e487affac531609a573c0faa57b9d823e1a0955 Author: Jérôme Hénin Date: Tue Apr 7 22:30:59 2020 +0200 Merge ABF output improvements from master commit 10588511133716b867e6428509e3bbb15ffbf918 Author: Jérôme Hénin Date: Mon Sep 16 19:12:17 2019 +0200 Merge name changes from master commit 50aaa9411c73b0e053f1b6a404d98b4bb5944072 Author: Jérôme Hénin Date: Wed Aug 1 19:14:32 2018 +0200 Reorder ABF function commit bd1723b5d9863eeab97ea79d424c0f3f3791554a Author: Jérôme Hénin Date: Wed Aug 1 19:01:41 2018 +0200 Merge improvements from linabf branch commit 687d6abf9681c5baaedcb5dfd05a37ce3e5516e2 Author: Jérôme Hénin Date: Fri Mar 9 11:08:46 2018 +0100 First impl. of lin-ABF commit 21628ab1656da9a084c2fa19d4547b6eef15fe9e Author: Jérôme Hénin Date: Fri Mar 2 17:02:40 2018 +0100 First version of Spring ABF (ABF-UI) commit cc9051e336266a6651c798a3d099e3776e1c6dc7 Author: Jérôme Hénin Date: Fri Mar 2 16:54:47 2018 +0100 Another layer of ABF-AR commit 6fea3593a6af316f83cde982d291466c6f59ed4c Author: Jérôme Hénin Date: Wed Feb 28 11:54:48 2018 +0100 WIP on ABF-AR: fullSamples commit 8f78f8f00ed59c22ca20c4f7159a58b568e571a7 Author: Jérôme Hénin Date: Mon Feb 26 17:45:06 2018 +0100 Another layer of WIP on ABF-AR commit 879a2cf2e9eadf67be308ddac5324d004a12e28a Author: Jérôme Hénin Date: Sun Feb 25 00:53:07 2018 +0100 WIP on ABF-AR commit 31853a6150e418bf1e7c1512365217d8e6f00655 Author: Jérôme Hénin Date: Fri Feb 23 16:49:29 2018 +0100 First version of ABF-AR for now, only 1D, not calculating expectation of biasing force commit 737741240673c42eaa36603d101d03eec54287a9 Author: Jérôme Hénin Date: Wed Aug 26 15:50:26 2020 +0200 as_array method for colvars_grid.py commit ac2c347c22513b97ddb573b61c82117295b76b11 Author: Jérôme Hénin Date: Tue Aug 11 18:03:15 2020 +0200 Add two methods ro colvars-grid python class scalar_array numerical_gradient (in 2d only) commit f4a3361c2e5172f2ed9a56c51021a49314b17f73 Author: Jérôme Hénin Date: Mon Aug 10 14:26:59 2020 +0200 Undo regression in grid::wrap() commit 2ba510e23471ba91fc4129f73c7aec6c1a379b93 Author: Jérôme Hénin Date: Mon Aug 10 14:26:39 2020 +0200 Do not create unnecessary temp grid when no CZAR commit 3b51d4e29939225b1d862ccde051dc1a0fe66160 Author: Jérôme Hénin Date: Wed Jul 8 16:45:46 2020 +0200 Fix grid_RMSD: special version for gradients commit b05294913c2dea2405ef31a402a1b36e3dd4d092 Author: Jérôme Hénin Date: Tue Jul 7 20:15:03 2020 +0200 Support for CZAR files in ABF_dataset class for post-analysis in Python commit 2550358286c4516753e4635d047d99514b9f3269 Author: Jérôme Hénin Date: Tue Jul 7 20:11:24 2020 +0200 Radius sampling for mwABF selection in d>1 Implemented at C++ level, exposing local_sample_count scripting command and deprecating bincount Results in simplified mwABF_selection script commit 748b0af2da3336bfcf0fe6054e530e39540be935 Author: Jérôme Hénin Date: Mon Jul 6 15:29:34 2020 +0200 Remove pointer aliasing when deleting ABF bias commit 9143a7f0b7926165e2f0b0141250b5eb10a91f3b Author: Jérôme Hénin Date: Sun Jul 5 23:23:43 2020 +0200 colvars_grid: add marginal_count function commit 9a9c5884ee75d2237ff3ff082772b478832018fe Author: Jérôme Hénin Date: Sun Jul 5 23:21:51 2020 +0200 Dashboard: warn about bogus PBC lengths from NAMD runs without PBC NAMD writes {1 1 1} as lengths into DCD files by default Offer to overwrite to avoid computing wrong colvar values commit 44bfdcd28337d8a8cf330b5ac9222bc88e0a52a2 Author: Jérôme Hénin Date: Wed Jul 1 20:59:15 2020 +0200 Print rmsd from global grad and pmf in mwABF commit 014c0d187053243a92b6cd16cf5ddc11814884a8 Author: Jérôme Hénin Date: Tue Jun 30 22:20:50 2020 +0200 Fix logic to maintain local CZAR data separately in shared eABF commit 40dc0cafe1ee33caf50397d9a2fa13a622018422 Author: Jérôme Hénin Date: Tue Jun 30 11:04:23 2020 +0200 Correctly integrate local PMF in shared ABF commit b8227f14ce34a00b9ca18433535498fd6f21aea0 Author: Jérôme Hénin Date: Tue Jun 30 11:03:28 2020 +0200 First version of ABF_dataset python class commit f295673f8e7db325510015c2b224c227a340eff2 Author: Jérôme Hénin Date: Fri Jun 26 17:06:53 2020 +0200 Correct logic for shared eABF/CZAR esp. for allocating the right grids at the right time commit c51c04c67becbd9e5cf964cea0fbcf14474d122a Author: Jérôme Hénin Date: Fri Jun 26 17:06:05 2020 +0200 Integrate can_accumulate_data into ABF commit 960dc341fee71ef89d40d6b4e32533e1502503ec Author: Jérôme Hénin Date: Fri Jun 26 14:05:36 2020 +0200 Allocate shared ABF grids at init time if possible avoids rare cases where shared_freq is more than output_freq and output is attempted from unallocated grids, leading to segfault. commit eaa8f70ac4c6c86be38fc03e86264b5feb2822e1 Author: Giacomo Fiorin Date: Thu Jun 25 15:31:44 2020 -0400 Add test for multiple run commands commit de7ab2bcff06084ca3e95cc83ac28227b0f6ba1e Author: Jérôme Hénin Date: Fri Jun 26 09:58:30 2020 +0200 mwABF script doc update commit 314a1aa78b1673ac1a69c07ff14b5a5cc580c3eb Author: Jérôme Hénin Date: Fri Jun 26 00:27:11 2020 +0200 Save local data per replica in shared ABF commit 69f53fd047c5c98869349c30c206d9946b39fed1 Author: Jérôme Hénin Date: Thu Jun 25 19:17:35 2020 +0200 Do not write ABF history files twice for same ts commit a3884dc8e84467b300ad9dfcfc11b31d101a28ec Author: Jérôme Hénin Date: Thu Jun 25 18:21:15 2020 +0200 Do not do sharing from mwABF selection script instead, expect shared ABF to be setup in Colvars config. This allows for better timing of ABF sharing vs. file IO (share first, write output files second). This is necessary now that output files are not written a second time at the end of a run if it coincides with outputFreq. The Tcl script is now only in charge of selection. It did not take any particular decision regarding sharing anyway, just duplicates the shareFreq behavior. commit 1c8f38db9ec7c5bfb67491850d6e6d0bee2bcdb3 Author: Jérôme Hénin Date: Wed Jun 24 15:11:43 2020 +0200 Fixes for colvars_grid.py commit 0449812363066f12d4bfa450aca3cb16a4a89bec Author: Jérôme Hénin Date: Thu Jun 18 22:45:24 2020 +0200 Improvements to mwABF - do sharing before force update for better mw-pABF - WIP on collecting and outputting local data per replica for variance analysis commit a08fe5d535538bf374940499235a4d296e42beef Author: Jérôme Hénin Date: Thu Jun 18 22:45:00 2020 +0200 Improvements to colvars_grid python class commit 5dc48fb32939ee43596179d4016380998f6ed607 Author: Jérôme Hénin Date: Wed Jun 17 13:32:40 2020 +0200 Write output files after a 0-step simulation commit 8f600c859a2471758a8e526de19fca6987cf003a Author: Jérôme Hénin Date: Thu May 28 23:56:37 2020 +0200 mwABF improvements - refactor mwABF selection scripts into single file - provide proc run_mwABF_selection with as many reasonable defaults as possible - toggle shareOn when sharing via scripting commands - make selection work in dimension > 1 - fix divide by zero bug when computing weights for selection - include "replicaUniformPatchGrids on" in selection script - share data for CZAR estimator if doing eABF --- colvartools/colvars_grid.py | 2 +- colvartools/mwABF_selection.tcl | 207 +++++ doc/colvars-refman-main.tex | 55 +- doc/cvscript-fix-modify.tex | 11 +- doc/cvscript-tcl.tex | 11 +- .../AutoDiff/test.colvars.out | 130 ++- .../AutoDiff/test.colvars.state.stripped | 18 +- .../AutoDiff/test.colvars.traj | 166 ++-- .../000_stringstates/AutoDiff/test.pmf | 84 +- .../000_stringstates/namd-version.txt | 6 +- .../005_Tcl_script/AutoDiff/test.Tcl.out | 14 +- .../005_Tcl_script/AutoDiff/test.colvars.out | 602 +++++++------ .../005_Tcl_script/AutoDiff/test.colvars.traj | 68 +- .../library/005_Tcl_script/AutoDiff/test.grad | 24 +- .../library/005_Tcl_script/AutoDiff/test.pmf | 24 +- .../abf_12-15_less.colvars | 1 - .../multiple_walker_abf/mwABF_selection.tcl | 1 + .../multiple_walker_abf/mw_independent.namd | 2 +- .../multiple_walker_abf/mw_selection.namd | 154 +--- .../multiple_walker_abf/run_mw_tests.sh | 15 +- .../sabf_12-15_less.colvars | 1 - .../multiple_walker_abf/selectionRules.tcl | 98 --- src/colvar.cpp | 27 +- src/colvaratoms.cpp | 2 +- src/colvarbias.cpp | 14 +- src/colvarbias.h | 3 +- src/colvarbias_abf.cpp | 801 +++++++++++------- src/colvarbias_abf.h | 114 ++- src/colvarcomp.cpp | 2 +- src/colvardeps.cpp | 5 + src/colvargrid.cpp | 121 ++- src/colvargrid.h | 456 ++++++++-- src/colvargrid_def.h | 3 + src/colvarscript_commands_bias.h | 23 +- vmd/cv_dashboard/cv_dashboard.tcl | 16 + 35 files changed, 2052 insertions(+), 1229 deletions(-) create mode 100644 colvartools/mwABF_selection.tcl create mode 120000 namd/tests/library/multiple_walker_abf/mwABF_selection.tcl delete mode 100644 namd/tests/library/multiple_walker_abf/selectionRules.tcl diff --git a/colvartools/colvars_grid.py b/colvartools/colvars_grid.py index 69d0b9097..4bc6609de 100644 --- a/colvartools/colvars_grid.py +++ b/colvartools/colvars_grid.py @@ -10,7 +10,7 @@ class colvars_grid: Gradient files contain as many data series as there are dimensions. Data can be obtained in array shape with as_array() - # Example: 3d surface plot for a 2d free energy surface ("PMF"), with contour plot at z=0 + # Example: 3d surface plot for a 2d free energy surface (PMF), with contour plot at z=0 from mpl_toolkits.mplot3d import Axes3D pmf = colvars_grid('run.pmf') fig = plt.figure() diff --git a/colvartools/mwABF_selection.tcl b/colvartools/mwABF_selection.tcl new file mode 100644 index 000000000..dd4b8d7db --- /dev/null +++ b/colvartools/mwABF_selection.tcl @@ -0,0 +1,207 @@ +# Scripts for multiple-walker ABF in NAMD with selection rules + +# NOTE: this script does not perform shared ABF - sharing must be enabled +# separately in the ABF configuration (shared on, sharedFreq xxx) + +# Source in NAMD config, and call: +# run_mwABF $cycleNum $selectFreq ($sampleRad $percentStop $biasName) + +# cycleNum Number of cycles of sharing ABF data to be performed +# selectFreq Number of timesteps between between rounds of replica selection +# NOTE: Total simulation time is cycleNum * selectFreq +# sampleRad Half-length (in bins) of the interval around the current position of a replica +# in which to count samples for selection purposes (default: 0) +# percentStop Relative difference (in percent) between lowest and highest sample counts +# below which no selection is performed (recommended range: 20 to 90, default: 40) +# biasName Name of the ABF bias to perform selection on (default: "abf1", Colvars default) + +# Author: Jeff Comer +# Comer et al. JCTC 2014 https://doi.org/10.1021/ct500874p +# With improvements by Jérôme Hénin + +# The following is needed for exchanging coordinates between NAMD replicas +# but cannot be "set from script", ie after a scripting command has been run +# So this file should be sourced before any script command (cv, run etc.) +replicaUniformPatchGrids on + + +proc run_mwABF_selection { cycleNum selectFreq {sampleRad 0} {percentStop 40} {biasName "abf1"}} { + + set repNum [numReplicas] + if { $repNum < 2 } { + print "Error: cannot perform shared ABF with fewer than 2 replicas.\nUse the +replicas command-line flag of NAMD (MPI build).\n\n" + exit 1 + } + set rep [myReplica] + set ncolvars 1 + # Parse ABF config to find out number of colvars + set config [cv bias $biasName getconfig] + if { [regexp -nocase {^\s*colvars\s+([^\s{}#]+)} $config match colvars] } { + set ncolvars [llength $colvars] + print "Shared ABF $biasName running on $ncolvars colvars: $colvars" + } + + set s "REPLICA_SETUP" + foreach var {cycleNum selectFreq repNum rep sampleRad percentStop} { + set s "$s $var [set $var]" + } + print "$s" + + for {set i 0} {$i < $cycleNum} {incr i} { + print "CYCLE $i" + # Run the steps until next selection + # Bias sharing is controlled by the ABF bias at the C++ level + run $selectFreq + + set count [cv bias $biasName local_sample_count $sampleRad] + + replicaBarrier + if {$rep > 0} { + ## Send the count to Replica 0. + replicaSend $count 0 + ## Receive the source and destination replicas. + set srcDestList [replicaRecv 0] + } else { + ## Build the weight list. + # The weight is the inverse count + 1 to avoid divide by 0 + set weightList [list [expr {1.0/($count + 1)}]] + print "REPLICA_COUNT $rep $count" + set countMin $count + set countMax $count + for {set r 1} {$r < $repNum} {incr r} { + set repCount [replicaRecv $r] + print "REPLICA_COUNT $r $repCount" + if {$repCount < $countMin} { + set countMin $repCount + } + if {$repCount > $countMax} { + set countMax $repCount + } + # The weight is the inverse count + 1 to avoid divide by 0 + lappend weightList [expr {1.0/($repCount + 1)}] + } + + ## Normalize the weight list. + set weightList [normalizeWeights $weightList] + set s "REPLICA_WEIGHT_LIST" + foreach w $weightList { + set s "$s [format " %.3g" $w]" + } + print $s + + print "REPLICA_MINMAX $countMin $countMax" + if {$countMin < 1} {set countMin 1} + set percentDif [expr {(100*($countMax - $countMin))/$countMax}] + + ## Generate the list of exchanges "srcDestList" + if {$percentDif < $percentStop} { + print "REPLICA_SELECTION_DISABLED $percentDif" + ## If the relative difference between the min and max counts + ##is less than the threshold, we don't do exchanges. + set srcDestList {} + } else { + print "REPLICA_SELECTION_ENABLED $percentDif" + set cloneList [resampleWalkers $weightList] + print "REPLICA_CLONE_LIST $cloneList" + + set srcDestList [resampleExchanges $cloneList] + } + ## Replica 0 sends the srcDestList to all other replicas, + ## so they know who to receive from and who to send to. + for {set r 1} {$r < $repNum} {incr r} { + replicaSend $srcDestList $r + } + } ;# End Replica 0 work + + ## Everyone should have an identical copy of the srcDestList. + if {[llength $srcDestList] > 0} { + print "REPLICA_SRC_DEST_LIST $srcDestList" + print "REPLICA_BARRIER" + replicaBarrier + + # Do the coordinate exchanges if this replica is the source or destination. + foreach srcDest $srcDestList { + set src [lindex $srcDest 0] + set dest [lindex $srcDest 1] + + if {$src == $rep} { + print "REPLICA_ATOM_SEND $dest" + replicaAtomSend $dest + } elseif {$dest == $rep} { + print "REPLICA_ATOM_RECV $src" + replicaAtomRecv $src + } + } + } else { + print "REPLICA_SRC_DEST_EMPTY" + } + } + return $i ;# Number of cycles completed +} + +proc normalizeWeights {weightList} { + ## Normalize the weights. + set weightSum 0.0 + foreach weight $weightList { + set weightSum [expr {$weightSum + $weight}] + } + + set r 0 + set wnList {} + foreach weight $weightList { + lappend wnList [expr {double($weight)/$weightSum}] + incr r + } + return $wnList +} + + +# Determine the number of clones for each walker. +# Weights must be normalized. +proc resampleWalkers {weightList} { + set num [llength $weightList] + + ## Get the number of clones for each walker. + set wbar(0) [lindex $weightList 0] + set u [expr {rand()}] + set cloneList [list [expr { int($num*$wbar(0)+$u) }] ] + for {set r 1} {$r < $num} {incr r} { + set r0 [expr {$r-1}] + set wbar($r) [expr {$wbar($r0) + [lindex $weightList $r]}] + lappend cloneList [expr {int($num*$wbar($r)+$u) - int($num*$wbar($r0)+$u) }] + } + return $cloneList +} + +# Determine the minimal exchanges that must be made to resample. +proc resampleExchanges {cloneList} { + ## Make a list of exchanges. + set cloneZeroList {} + set cloneMultList {} + set r 0 + foreach cloneNum $cloneList { + if {$cloneNum == 0} { lappend cloneZeroList $r } + if {$cloneNum > 1} { lappend cloneMultList $r } + incr r + } + + # Is nothing cloned? + if {[llength $cloneZeroList] == 0} { + return {} + } + + ## Walkers cloned multiple times are copied to walkers cloned zero times. + ## Make the list of exchanges srcDestList + set srcDestList {} + set zeroInd 0 + foreach mult $cloneMultList { + # We get one clone just by leaving the walker where it is. + set extraNum [expr {[lindex $cloneList $mult]-1}] + for {set j 0} {$j < $extraNum} {incr j} { + set dest [lindex $cloneZeroList $zeroInd] + lappend srcDestList [list $mult $dest] + incr zeroInd + } + } + return $srcDestList +} diff --git a/doc/colvars-refman-main.tex b/doc/colvars-refman-main.tex index b0946fd02..c57b67d7b 100644 --- a/doc/colvars-refman-main.tex +++ b/doc/colvars-refman-main.tex @@ -5552,35 +5552,52 @@ \end{itemize} \cvnamdonly{ -\cvsubsubsec{Multiple-replica ABF}{sec:colvarbias_abf_shared} +\cvsubsubsec{Multiple-walker ABF}{sec:colvarbias_abf_shared} \label{sec:mw-ABF} +This implements the multiple-walker ABF scheme described in \cite{Minoukadeh2010}. The reference for this +implementation is \cite{Comer2014c}. +This feature requires that \MDENGINE be compiled and executed with multiple-replica +support. + +If \refkey{shared}{abf|shared} is enabled, the total force samples will be synchronized among all replicas +at intervals defined by \refkey{sharedFreq}{abf|sharedFreq}. +Each replica maintains a separate buffer of total force samples that determine the biasing force. +Every \refkey{sharedFreq}{abf|sharedFreq} steps, the replicas communicate the samples that have been gathered since the last synchronization time, ensuring all replicas apply a similar biasing force. +Thus, it is as if total force samples among all replicas are gathered in a single shared buffer. +Shared ABF allows all replicas to benefit from the sampling done by other replicas and can lead to faster convergence of the biasing force. + +An implementation of the selection mechanism described in \cite{Comer2014c} is provided as a set of Tcl scripts in \texttt{colvartools/mwABF\_selection.tcl} in the Colvars repository. +Compared with the initial implementation, the current implementation supports selection on sets of up to 3 colvars. +To reduce noise, samples can be counted in a hybercube around the current bin (see Tcl script for details). +A set of example NAMD inputs can be found in the \texttt{namd/tests/library/multiple\_walker\_abf} directory. + + +\paragraph{Output files of multiple-walker ABF.} +In multiple-walker ABF runs, each walker now outputs gradient and count files containing only data collected locally (this is a change in version 2024-01-06). +In addition, the first walker outputs the collected data using the common prefix, with an additional \texttt{".all"} string (i.e. file names ending with \texttt{".all.count"}, \texttt{".all.grad"} etc.). + + \begin{itemize} -\item \keydef{shared}{\texttt{abf}}{% - Apply multiple-replica ABF, sharing force samples among the replicas?} +\item + \labelkey{abf|shared} + \keydef{shared}{\texttt{abf}}{% + Apply multiple-walker ABF, sharing force samples among the replicas?} {boolean} {\texttt{no}} - { This is command requires that NAMD be compiled and executed with multiple-replica - support. - If \texttt{shared} is set to yes, the total force samples will be synchronized among all replicas - at intervals defined by \texttt{sharedFreq}. - This implements the multiple-walker ABF scheme described in \cite{Minoukadeh2010}; this - implementation is documented in \cite{Comer2014c}. - Thus, it is as if total force samples among all replicas are - gathered in a single shared buffer, which why the algorithm is referred to as shared ABF. - Shared ABF allows all replicas to benefit from the sampling done by other replicas and can lead to faster convergence of the biasing force. + { + Enable sharing of ABF data between replicas, in a multi-replica simulation. } -\item \keydef{sharedFreq}{\texttt{abf}}{% +\item + \labelkey{abf|sharedFreq} + \keydef{sharedFreq}{\texttt{abf}}{% Frequency (in timesteps) at which force samples are synchronized among the replicas} - {positive integer} + {positive integer or zero} {\refkey{outputFreq}{colvarbias|outputFreq}} { - In the current implementation of shared ABF, each replica maintains a separate - buffer of total force samples that determine the biasing force. - Every \texttt{sharedFreq} steps, the replicas communicate the samples that - have been gathered since the last synchronization time, ensuring all replicas - apply a similar biasing force. + Every \texttt{sharedFreq} steps, the replicas communicate the samples that have been gathered since the last synchronization time. + Because this forces the replicas to synchronize, small values of this parameter may slightly reduce the overall throughput, especially on systems with unequal performance. } \end{itemize} } diff --git a/doc/cvscript-fix-modify.tex b/doc/cvscript-fix-modify.tex index e7fab8a7d..c1ab5dcde 100644 --- a/doc/cvscript-fix-modify.tex +++ b/doc/cvscript-fix-modify.tex @@ -595,7 +595,7 @@ \begin{mdexampleinput}{} \texttt{\textbf{fix\_modify Colvars bias name bin}} \\ -\-~~~~\texttt{Get the current grid bin index (1D ABF only for now)} +\-~~~~\texttt{Get the current grid bin index (flattened if more than 1d)} \\ \-~~~~\texttt{\textbf{Returns}} \\ @@ -615,6 +615,15 @@ \-~~~~\texttt{samples : integer - Number of samples} \end{mdexampleinput} \begin{mdexampleinput}{} +\texttt{\textbf{fix\_modify Colvars bias name local\_sample\_count [radius]}} +\\ +\-~~~~\texttt{Get the number of samples around the current binsamples : integer - Number of samples} +\\ +\-~~~~\texttt{\textbf{Parameters}} +\\ +\-~~~~\texttt{radius : integer - Sum over radius bins around current bin (optional)} +\end{mdexampleinput} +\begin{mdexampleinput}{} \texttt{\textbf{fix\_modify Colvars bias name binnum}} \\ \-~~~~\texttt{Get the total number of grid points of this bias (1D ABF only for now)} diff --git a/doc/cvscript-tcl.tex b/doc/cvscript-tcl.tex index 15f23e271..1092a030c 100644 --- a/doc/cvscript-tcl.tex +++ b/doc/cvscript-tcl.tex @@ -595,7 +595,7 @@ \begin{mdexampleinput}{} \texttt{\textbf{cv bias name bin}} \\ -\-~~~~\texttt{Get the current grid bin index (1D ABF only for now)} +\-~~~~\texttt{Get the current grid bin index (flattened if more than 1d)} \\ \-~~~~\texttt{\textbf{Returns}} \\ @@ -615,6 +615,15 @@ \-~~~~\texttt{samples : integer - Number of samples} \end{mdexampleinput} \begin{mdexampleinput}{} +\texttt{\textbf{cv bias name local\_sample\_count [radius]}} +\\ +\-~~~~\texttt{Get the number of samples around the current binsamples : integer - Number of samples} +\\ +\-~~~~\texttt{\textbf{Parameters}} +\\ +\-~~~~\texttt{radius : integer - Sum over radius bins around current bin (optional)} +\end{mdexampleinput} +\begin{mdexampleinput}{} \texttt{\textbf{cv bias name binnum}} \\ \-~~~~\texttt{Get the total number of grid points of this bias (1D ABF only for now)} diff --git a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.out b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.out index bf8e1ba9b..da43847c9 100644 --- a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.out +++ b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.out @@ -1,13 +1,12 @@ colvars: ---------------------------------------------------------------------- colvars: Please cite Fiorin et al, Mol Phys 2013: -colvars: https://dx.doi.org/10.1080/00268976.2013.813594 -colvars: in any publication based on this calculation. +colvars: https://doi.org/10.1080/00268976.2013.813594 +colvars: as well as all other papers listed below for individual features used. colvars: SMP parallelism is enabled; if needed, use "smp off" to override this. colvars: This version was built with the C++11 standard or higher. +colvars: Redefining the Tcl "cv" command to the new script interface. colvars: The restart output state file will be "test.tmp.colvars.state". colvars: The final output state file will be "test.colvars.state". -colvars: Opening trajectory file "test.colvars.traj". -colvars: Redefining the Tcl "cv" command to the new script interface. colvars: ---------------------------------------------------------------------- colvars: Reading new configuration from file "test.in": colvars: # units = "" [default] @@ -46,6 +45,7 @@ colvars: # colvarsTrajFrequency = 1 colvars: # colvarsRestartFrequency = 10 colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "one" @@ -59,28 +59,28 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "main". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group3" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "main". colvars: # printAtomIDs = off [default] colvars: Atom group "main" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. colvars: Initializing atom group "ref". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group4" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "ref". colvars: # printAtomIDs = off [default] @@ -110,18 +110,20 @@ colvars: ---------------------------------------------------------------------- colvars: Initializing a new "abf" instance. colvars: # name = "abf1" [default] colvars: # colvars = { one } +colvars: # stepZeroData = off [default] colvars: # outputEnergy = off [default] colvars: # outputFreq = 10 [default] colvars: # timeStepFactor = 1 [default] colvars: WARNING: ABF should not be run without a thermostat or at 0 Kelvin! colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 1 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -133,14 +135,35 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 8 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: +colvars: SUMMARY OF COLVARS FEATURES USED SO FAR AND THEIR CITATIONS: +colvars: +colvars: - Colvars module: +colvars: - Colvars-NAMD interface: +colvars: - Optimal rotation via flexible fitting: +colvars: - distanceZ colvar component: +colvars: Fiorin2013 https://doi.org/10.1080/00268976.2013.813594 +colvars: +colvars: - ABF colvar bias implementation: +colvars: - Internal-forces free energy estimator: +colvars: Henin2010 https://doi.org/10.1021/ct9004432 +colvars: +colvars: - NAMD engine: +colvars: - Scalable center-of-mass computation (NAMD): +colvars: Phillips2020 https://doi.org/10.1063/5.0014475 +colvars: +colvars: updating target temperature (T = 0 K). colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 8 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 0 K). +colvars: Current simulation parameters: initial step = 0, integration timestep = 1 +colvars: Updating atomic parameters (masses, charges, etc). colvars: Re-initialized atom group for variable "one":0/0. 4 atoms: total mass = 54.028, total charge = -0.4. colvars: Re-initialized atom group for variable "one":0/1. 4 atoms: total mass = 54.028, total charge = -0.4. colvars: The restart output state file will be "test.tmp.colvars.state". -colvars: The final output state file will be "test.colvars.state". -colvars: Prepared sample and gradient buffers at step 0. colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Saving collective variables state to "test.tmp.colvars.state". @@ -187,6 +210,7 @@ colvars: # colvarsTrajFrequency = 1 colvars: # colvarsRestartFrequency = 10 colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "one" @@ -200,32 +224,32 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "main". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group3" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "main". colvars: # printAtomIDs = off [default] -colvars: Atom group "main" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "main" defined with 4 atoms requested. colvars: Initializing atom group "ref". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group4" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "ref". colvars: # printAtomIDs = off [default] -colvars: Atom group "ref" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "ref" defined with 4 atoms requested. colvars: # axis = ( 0 , 0 , 1 ) [default] colvars: # oneSiteSystemForce = off [default] colvars: # oneSiteTotalForce = off [default] @@ -251,18 +275,20 @@ colvars: ---------------------------------------------------------------------- colvars: Initializing a new "abf" instance. colvars: # name = "abf1" [default] colvars: # colvars = { one } +colvars: # stepZeroData = off [default] colvars: # outputEnergy = off [default] colvars: # outputFreq = 10 [default] colvars: # timeStepFactor = 1 [default] colvars: WARNING: ABF should not be run without a thermostat or at 0 Kelvin! colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 1 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -274,7 +300,9 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 8 atoms in total). -colvars: Prepared sample and gradient buffers at step 20. +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 0 K). +colvars: The restart output state file will be "test.tmp.colvars.state". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Saving collective variables state to "test.tmp.colvars.state". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". @@ -321,6 +349,7 @@ colvars: # colvarsTrajFrequency = 1 colvars: # colvarsRestartFrequency = 10 colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "one" @@ -334,32 +363,32 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "main". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group3" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "main". colvars: # printAtomIDs = off [default] -colvars: Atom group "main" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "main" defined with 4 atoms requested. colvars: Initializing atom group "ref". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group4" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "ref". colvars: # printAtomIDs = off [default] -colvars: Atom group "ref" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "ref" defined with 4 atoms requested. colvars: # axis = ( 0 , 0 , 1 ) [default] colvars: # oneSiteSystemForce = off [default] colvars: # oneSiteTotalForce = off [default] @@ -385,18 +414,20 @@ colvars: ---------------------------------------------------------------------- colvars: Initializing a new "abf" instance. colvars: # name = "abf1" [default] colvars: # colvars = { one } +colvars: # stepZeroData = off [default] colvars: # outputEnergy = off [default] colvars: # outputFreq = 10 [default] colvars: # timeStepFactor = 1 [default] colvars: WARNING: ABF should not be run without a thermostat or at 0 Kelvin! colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 1 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -408,12 +439,14 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 8 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 0 K). colvars: ---------------------------------------------------------------------- -colvars: Loading state from input buffer. +colvars: Loading state from formatted string. colvars: Restarting collective variable "one" from value: 2.12037 -colvars: Restarting abf bias "abf1" from step number 20. +colvars: Restarted abf bias "abf1" with step number 20. colvars: ---------------------------------------------------------------------- -colvars: Prepared sample and gradient buffers at step 20. +colvars: The restart output state file will be "test.tmp.colvars.state". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Saving collective variables state to "test.tmp.colvars.state". @@ -459,6 +492,7 @@ colvars: # colvarsTrajFrequency = 1 colvars: # colvarsRestartFrequency = 10 colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "one" @@ -472,32 +506,32 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "main". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group3" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "main". colvars: # printAtomIDs = off [default] -colvars: Atom group "main" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "main" defined with 4 atoms requested. colvars: Initializing atom group "ref". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "group4" colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "ref". colvars: # printAtomIDs = off [default] -colvars: Atom group "ref" defined with 4 atoms requested: total mass = 54.028, total charge = -0.4. +colvars: Atom group "ref" defined with 4 atoms requested. colvars: # axis = ( 0 , 0 , 1 ) [default] colvars: # oneSiteSystemForce = off [default] colvars: # oneSiteTotalForce = off [default] @@ -523,18 +557,20 @@ colvars: ---------------------------------------------------------------------- colvars: Initializing a new "abf" instance. colvars: # name = "abf1" [default] colvars: # colvars = { one } +colvars: # stepZeroData = off [default] colvars: # outputEnergy = off [default] colvars: # outputFreq = 10 [default] colvars: # timeStepFactor = 1 [default] colvars: WARNING: ABF should not be run without a thermostat or at 0 Kelvin! colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 1 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -546,8 +582,10 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 8 atoms in total). -colvars: Restarting abf bias "abf1" from step number 40. -colvars: Prepared sample and gradient buffers at step 40. +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 0 K). +colvars: Restarted abf bias "abf1" with step number 40. +colvars: The restart output state file will be "test.tmp.colvars.state". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". colvars: Saving collective variables state to "test.tmp.colvars.state". colvars: Synchronizing (emptying the buffer of) trajectory file "test.colvars.traj". diff --git a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.state.stripped b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.state.stripped index 9cdeb4354..902e2f470 100644 --- a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.state.stripped +++ b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.state.stripped @@ -5,27 +5,27 @@ configuration { colvar { name one - x 1.93953985558713e+00 + x 1.9395398554137 } abf { configuration { - step 60 - name abf1 +step 60 +name abf1 } samples - 1 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 11 49 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 gradient - 3.20950796166942e+01 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 - 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 - 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 -1.54059351654355e+01 - -1.77818813718884e+01 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 - 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 0.00000000000000e+00 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 -15.40593524979 + -17.781881376907 0 0 0 0 0 0 0 + 0 0 0 0 0 0 0 0 } diff --git a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.traj b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.traj index 3bd3feed3..a65642d66 100644 --- a/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.traj +++ b/namd/tests/interface/000_stringstates/AutoDiff/test.colvars.traj @@ -1,88 +1,88 @@ # step one fa_one 0 2.13236862431278e+00 0.00000000000000e+00 - 1 2.13245751463817e+00 -2.89496599117013e+01 - 2 2.13267749923815e+00 -3.31817108699520e+01 - 3 2.13306271104622e+00 -3.67376757993221e+01 - 4 2.13362553120203e+00 -3.96044965685618e+01 - 5 2.13435489235458e+00 -4.17550404792339e+01 - 6 2.13521727275388e+00 -4.31864384367706e+01 - 7 2.13615933633749e+00 -4.39214019856548e+01 - 8 2.13711157280064e+00 -4.40035034823288e+01 - 9 2.13799240955077e+00 -4.34912898697864e+01 - 10 2.13871247289936e+00 -4.24532563931769e+01 - 11 2.13917889857984e+00 -4.09645185575262e+01 - 12 2.13929973498283e+00 -3.91051042320417e+01 - 13 2.13898846494078e+00 -3.69590309335845e+01 - 14 2.13816852528690e+00 -3.46131556629380e+01 - 15 2.13677750977939e+00 -3.21550548077171e+01 - 16 2.13477057601376e+00 -2.96696321389358e+01 - 17 2.13212255550928e+00 -2.72348223934214e+01 - 18 2.12882840547568e+00 -2.49171916472515e+01 - 19 2.12490189621160e+00 -2.27683351162044e+01 - 20 2.12037277644266e+00 -2.08229577360183e+01 + 1 2.13245751463822e+00 -2.89496599179517e+01 + 2 2.13267749923821e+00 -3.31817108750264e+01 + 3 2.13306271104623e+00 -3.67376758031010e+01 + 4 2.13362553120190e+00 -3.96044965709237e+01 + 5 2.13435489235419e+00 -4.17550404800491e+01 + 6 2.13521727275308e+00 -4.31864384359511e+01 + 7 2.13615933633610e+00 -4.39214019832102e+01 + 8 2.13711157279848e+00 -4.40035034783994e+01 + 9 2.13799240954768e+00 -4.34912898646446e+01 + 10 2.13871247289520e+00 -4.24532563871765e+01 + 11 2.13917889857452e+00 -4.09645185510267e+01 + 12 2.13929973497631e+00 -3.91051042253309e+01 + 13 2.13898846493305e+00 -3.69590309268313e+01 + 14 2.13816852527797e+00 -3.46131556561760e+01 + 15 2.13677750976923e+00 -3.21550548008751e+01 + 16 2.13477057600232e+00 -2.96696321318850e+01 + 17 2.13212255549648e+00 -2.72348223860342e+01 + 18 2.12882840546139e+00 -2.49171916394546e+01 + 19 2.12490189619572e+00 -2.27683351080122e+01 + 20 2.12037277642510e+00 -2.08229577275481e+01 # step one fa_one - 20 2.12037277644266e+00 0.00000000000000e+00 - 21 2.11544271623815e+00 1.55685198258614e+01 - 22 2.11053109498525e+00 1.49732643651787e+01 - 23 2.10564809931018e+00 1.43571942152212e+01 - 24 2.10079970229767e+00 1.37987285480790e+01 - 25 2.09598736284211e+00 1.33331173109700e+01 - 26 2.09120837141256e+00 1.29742647577597e+01 - 27 2.08645727544484e+00 1.27169689706410e+01 - 28 2.08172838292490e+00 1.25377732682915e+01 - 29 2.07701914059296e+00 1.23968050016051e+01 - 30 2.07233397074510e+00 1.22413810067053e+01 - 31 2.06768794096187e+00 1.20116466162023e+01 - 32 2.06310948669217e+00 1.16480621935819e+01 - 33 2.05864143962384e+00 1.10997205765136e+01 - 34 2.05433986586549e+00 1.03319813367191e+01 - 35 2.05027079226797e+00 9.33122313524745e+00 - 36 2.04650549391532e+00 8.10539508119449e+00 - 37 2.04311550284503e+00 6.68000008657994e+00 - 38 2.04016847954816e+00 5.09109927809035e+00 - 39 2.03772551925461e+00 3.37821358816879e+00 - 40 2.03583969535811e+00 1.57963799786840e+00 + 20 2.12037277642510e+00 0.00000000000000e+00 + 21 2.11544271621882e+00 1.55685198357011e+01 + 22 2.11053109496424e+00 1.49732643721063e+01 + 23 2.10564809928765e+00 1.43571942188257e+01 + 24 2.10079970227382e+00 1.37987285482642e+01 + 25 2.09598736281719e+00 1.33331173079539e+01 + 26 2.09120837138683e+00 1.29742647520364e+01 + 27 2.08645727541851e+00 1.27169689629116e+01 + 28 2.08172838289812e+00 1.25377732593621e+01 + 29 2.07701914056579e+00 1.23968049923061e+01 + 30 2.07233397071747e+00 1.22413809978373e+01 + 31 2.06768794093357e+00 1.20116466085310e+01 + 32 2.06310948666284e+00 1.16480621878565e+01 + 33 2.05864143959295e+00 1.10997205734801e+01 + 34 2.05433986583228e+00 1.03319813370927e+01 + 35 2.05027079223152e+00 9.33122313964807e+00 + 36 2.04650549387451e+00 8.10539509006982e+00 + 37 2.04311550279864e+00 6.68000010013434e+00 + 38 2.04016847949489e+00 5.09109929626650e+00 + 39 2.03772551919320e+00 3.37821361066347e+00 + 40 2.03583969528738e+00 1.57963802415380e+00 # step one fa_one - 20 2.03583969535811e+00 -2.08229577360183e+01 - 21 2.03438314391287e+00 -2.16143395489082e+01 - 22 2.03323634182338e+00 -2.25233835188036e+01 - 23 2.03244195151392e+00 -2.35126646834306e+01 - 24 2.03202876537072e+00 -2.45381833176116e+01 - 25 2.03200708220403e+00 -2.55492789140476e+01 - 26 2.03236463490451e+00 -2.64910236828290e+01 - 27 2.03306411080377e+00 -2.73086633692016e+01 - 28 2.03404310563579e+00 -2.79532040855993e+01 - 29 2.03521691814417e+00 -2.83869421335027e+01 - 30 2.03648407621760e+00 -2.85878373489444e+01 - 31 2.03773389344912e+00 -2.85517208499622e+01 - 32 2.03885493501486e+00 -2.82918958206758e+01 - 33 2.03974309290751e+00 -2.78362553815613e+01 - 34 2.04030808991177e+00 -2.72225823535986e+01 - 35 2.04047757093319e+00 -2.64929683631305e+01 - 36 2.04019848933884e+00 -2.56884837644217e+01 - 37 2.03943606554932e+00 -2.48450516261266e+01 - 38 2.03817105695231e+00 -2.39911203526063e+01 - 39 2.03639641404859e+00 -2.31474142163199e+01 - 40 2.03411446673546e+00 -2.23285614789767e+01 + 20 2.03583969528738e+00 -2.08229577275480e+01 + 21 2.03438314383169e+00 -2.16143395365368e+01 + 22 2.03323634173075e+00 -2.25233835035229e+01 + 23 2.03244195140922e+00 -2.35126646663931e+01 + 24 2.03202876525377e+00 -2.45381833000924e+01 + 25 2.03200708207517e+00 -2.55492788973770e+01 + 26 2.03236463476459e+00 -2.64910236683028e+01 + 27 2.03306411065418e+00 -2.73086633579748e+01 + 28 2.03404310547835e+00 -2.79532040785819e+01 + 29 2.03521691798102e+00 -2.83869421312746e+01 + 30 2.03648407605107e+00 -2.85878373517011e+01 + 31 2.03773389328150e+00 -2.85517208575191e+01 + 32 2.03885493484829e+00 -2.82918958325225e+01 + 33 2.03974309274381e+00 -2.78362553969507e+01 + 34 2.04030808975232e+00 -2.72225823716478e+01 + 35 2.04047757077894e+00 -2.64929683829132e+01 + 36 2.04019848919025e+00 -2.56884837850437e+01 + 37 2.03943606540642e+00 -2.48450516467789e+01 + 38 2.03817105681471e+00 -2.39911203725986e+01 + 39 2.03639641391554e+00 -2.31474142350967e+01 + 40 2.03411446660596e+00 -2.23285614961233e+01 # step one fa_one - 40 2.03411446673546e+00 -2.23285614789767e+01 - 41 2.03133542224395e+00 -2.15457546590573e+01 - 42 2.02807738961372e+00 -2.08094244338103e+01 - 43 2.02436747932636e+00 -2.01309248710146e+01 - 44 2.02024294497527e+00 -1.95225111682743e+01 - 45 2.01575128278491e+00 -1.89957642980454e+01 - 46 2.01094849876049e+00 -1.85590777458691e+01 - 47 2.00589538855945e+00 -1.82151988151808e+01 - 48 2.00065242663384e+00 -1.79598284670796e+01 - 49 1.99527441618993e+00 0.00000000000000e+00 - 50 1.99007868263627e+00 -1.17687075273328e+01 - 51 1.98491160977150e+00 -1.26939462247035e+01 - 52 1.97977362511436e+00 -1.33199729391028e+01 - 53 1.97466058715707e+00 -1.37226438617414e+01 - 54 1.96956698074140e+00 -1.39735633609260e+01 - 55 1.96448934176169e+00 -1.41453862797532e+01 - 56 1.95942870840908e+00 -1.43022388973467e+01 - 57 1.95439148285722e+00 -1.44911544810089e+01 - 58 1.94938864094287e+00 -1.47378188783809e+01 - 59 1.94443365555237e+00 -1.50468045274783e+01 - 60 1.93953985558713e+00 -1.54059351654355e+01 + 40 2.03411446660596e+00 -2.23285614961230e+01 + 41 2.03133542211679e+00 -2.15457546742956e+01 + 42 2.02807738948757e+00 -2.08094244469945e+01 + 43 2.02436747919983e+00 -2.01309248821246e+01 + 44 2.02024294484702e+00 -1.95225111774162e+01 + 45 2.01575128265373e+00 -1.89957643054531e+01 + 46 2.01094849862541e+00 -1.85590777519031e+01 + 47 2.00589538841980e+00 -1.82151988203128e+01 + 48 2.00065242648937e+00 -1.79598284718619e+01 + 49 1.99527441604081e+00 0.00000000000000e+00 + 50 1.99007868248320e+00 -1.17687075716583e+01 + 51 1.98491160961488e+00 -1.26939462823049e+01 + 52 1.97977362495473e+00 -1.33199730082051e+01 + 53 1.97466058699500e+00 -1.37226439402165e+01 + 54 1.96956698057744e+00 -1.39735634464329e+01 + 55 1.96448934159629e+00 -1.41453863699019e+01 + 56 1.95942870824248e+00 -1.43022389898552e+01 + 57 1.95439148268945e+00 -1.44911545738228e+01 + 58 1.94938864077374e+00 -1.47378189697373e+01 + 59 1.94443365538142e+00 -1.50468046159132e+01 + 60 1.93953985541366e+00 -1.54059352497904e+01 diff --git a/namd/tests/interface/000_stringstates/AutoDiff/test.pmf b/namd/tests/interface/000_stringstates/AutoDiff/test.pmf index 5099f058a..4e0fcd82e 100644 --- a/namd/tests/interface/000_stringstates/AutoDiff/test.pmf +++ b/namd/tests/interface/000_stringstates/AutoDiff/test.pmf @@ -1,44 +1,44 @@ # 1 -# -10.25 0.5 41 0 +# -1.02500000000000e+01 5.00000000000000e-01 41 0 - -10 0.546368 - -9.5 16.5939 - -9 16.5939 - -8.5 16.5939 - -8 16.5939 - -7.5 16.5939 - -7 16.5939 - -6.5 16.5939 - -6 16.5939 - -5.5 16.5939 - -5 16.5939 - -4.5 16.5939 - -4 16.5939 - -3.5 16.5939 - -3 16.5939 - -2.5 16.5939 - -2 16.5939 - -1.5 16.5939 - -1 16.5939 - -0.5 16.5939 - 0 16.5939 - 0.5 16.5939 - 1 16.5939 - 1.5 16.5939 - 2 8.89094 - 2.5 0 - 3 0 - 3.5 0 - 4 0 - 4.5 0 - 5 0 - 5.5 0 - 6 0 - 6.5 0 - 7 0 - 7.5 0 - 8 0 - 8.5 0 - 9 0 - 9.5 0 - 10 0 + -1.00000000000000e+01 1.65939083133488e+01 + -9.50000000000000e+00 1.65939083133488e+01 + -9.00000000000000e+00 1.65939083133488e+01 + -8.50000000000000e+00 1.65939083133488e+01 + -8.00000000000000e+00 1.65939083133488e+01 + -7.50000000000000e+00 1.65939083133488e+01 + -7.00000000000000e+00 1.65939083133488e+01 + -6.50000000000000e+00 1.65939083133488e+01 + -6.00000000000000e+00 1.65939083133488e+01 + -5.50000000000000e+00 1.65939083133488e+01 + -5.00000000000000e+00 1.65939083133488e+01 + -4.50000000000000e+00 1.65939083133488e+01 + -4.00000000000000e+00 1.65939083133488e+01 + -3.50000000000000e+00 1.65939083133488e+01 + -3.00000000000000e+00 1.65939083133488e+01 + -2.50000000000000e+00 1.65939083133488e+01 + -2.00000000000000e+00 1.65939083133488e+01 + -1.50000000000000e+00 1.65939083133488e+01 + -1.00000000000000e+00 1.65939083133488e+01 + -5.00000000000000e-01 1.65939083133488e+01 + 0.00000000000000e+00 1.65939083133488e+01 + 5.00000000000000e-01 1.65939083133488e+01 + 1.00000000000000e+00 1.65939083133488e+01 + 1.50000000000000e+00 1.65939083133488e+01 + 2.00000000000000e+00 8.89094068845353e+00 + 2.50000000000000e+00 0.00000000000000e+00 + 3.00000000000000e+00 0.00000000000000e+00 + 3.50000000000000e+00 0.00000000000000e+00 + 4.00000000000000e+00 0.00000000000000e+00 + 4.50000000000000e+00 0.00000000000000e+00 + 5.00000000000000e+00 0.00000000000000e+00 + 5.50000000000000e+00 0.00000000000000e+00 + 6.00000000000000e+00 0.00000000000000e+00 + 6.50000000000000e+00 0.00000000000000e+00 + 7.00000000000000e+00 0.00000000000000e+00 + 7.50000000000000e+00 0.00000000000000e+00 + 8.00000000000000e+00 0.00000000000000e+00 + 8.50000000000000e+00 0.00000000000000e+00 + 9.00000000000000e+00 0.00000000000000e+00 + 9.50000000000000e+00 0.00000000000000e+00 + 1.00000000000000e+01 0.00000000000000e+00 diff --git a/namd/tests/interface/000_stringstates/namd-version.txt b/namd/tests/interface/000_stringstates/namd-version.txt index e2b19c981..9cc9ad5be 100644 --- a/namd/tests/interface/000_stringstates/namd-version.txt +++ b/namd/tests/interface/000_stringstates/namd-version.txt @@ -1,3 +1,3 @@ -Info: NAMD 2.14b1 for Linux-x86_64-multicore -colvars: Initializing the collective variables module, version "2020-04-09". -colvars: Using NAMD interface, version "2020-04-07". +Info: NAMD 3.0b4 for Linux-x86_64-multicore +colvars: Initializing the collective variables module, version 2023-10-03. +colvars: Using NAMD interface, version "2023-10-03". diff --git a/namd/tests/library/005_Tcl_script/AutoDiff/test.Tcl.out b/namd/tests/library/005_Tcl_script/AutoDiff/test.Tcl.out index f8428a56e..5e71fe400 100644 --- a/namd/tests/library/005_Tcl_script/AutoDiff/test.Tcl.out +++ b/namd/tests/library/005_Tcl_script/AutoDiff/test.Tcl.out @@ -1,4 +1,4 @@ -TCL: Running Colvars version 2020-04-09 +TCL: Running Colvars version 2023-10-03 TCL: Running cv units real TCL: TCL: d = 0.00000000000000e+00 @@ -13,7 +13,7 @@ TCL: E(harmonic1) = 18000 TCL: Original numsteps 20 will be ignored. TCL: Running for 20 steps TCL: d = 1.28684277602119e+01 -TCL: E(abf1) = 0 +TCL: E(abf1) = 32.885 TCL: E(harmonic1) = 94.2708 TCL: Running cv bias harm2 getconfig TCL: @@ -25,12 +25,12 @@ TCL: TCL: Running cv bias harmonic1 delete TCL: TCL: Running cv bias harm2 energy -TCL: 17.8166 +TCL: 17.8044 TCL: Running cv bias harm2 update -TCL: 17.8166 +TCL: 1.78043835964990e+01 TCL: Running for 20 steps TCL: Running cv bias harm2 energy -TCL: 13.4817 +TCL: 13.4812 TCL: Running cv bias harm2 state TCL: TCL: Running cv bias harm2 delete @@ -42,8 +42,8 @@ TCL: Atom groups: {3 } {98 } TCL: Gradients computed? (0/1): 1 TCL: Atom IDs: TCL: Gradients: -TCL: Applied force: 2.59631860820309e+01 -TCL: Total force: -1.36757725642476e+02 +TCL: Applied force: 2.59626762337285e+01 +TCL: Total force: -1.36804624327918e+02 TCL: Re-saving the state with a prefix: TCL: Re-saving the state with an explicit name: TCL: And reload: diff --git a/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.out b/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.out index d3135664a..ddbfc19bb 100644 --- a/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.out +++ b/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.out @@ -1,9 +1,10 @@ colvars: ---------------------------------------------------------------------- colvars: Please cite Fiorin et al, Mol Phys 2013: -colvars: https://dx.doi.org/10.1080/00268976.2013.813594 -colvars: in any publication based on this calculation. +colvars: https://doi.org/10.1080/00268976.2013.813594 +colvars: as well as all other papers listed below for individual features used. colvars: SMP parallelism is enabled; if needed, use "smp off" to override this. colvars: This version was built with the C++11 standard or higher. +colvars: Redefining the Tcl "cv" command to the new script interface. colvars: ---------------------------------------------------------------------- colvars: Reading new configuration from file "test.in": colvars: # units = "" [default] @@ -12,6 +13,7 @@ colvars: # colvarsTrajFrequency = 2 colvars: # colvarsRestartFrequency = 0 [default] colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "d" @@ -25,28 +27,28 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "group1". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "" [default] colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "group1". colvars: # printAtomIDs = off [default] colvars: Atom group "group1" defined with 1 atoms requested: total mass = 12.011, total charge = -0.1. colvars: Initializing atom group "group2". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "" [default] colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "group2". colvars: # printAtomIDs = off [default] @@ -80,13 +82,14 @@ colvars: # outputEnergy = off [default] colvars: # outputFreq = 60 colvars: # timeStepFactor = 1 [default] colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 0 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -104,6 +107,7 @@ colvars: # centers = { 12 } colvars: # targetCenters = { 12 } [default] colvars: # outputCenters = off [default] colvars: # forceConstant = 10 +colvars: # decoupling = off [default] colvars: # targetForceConstant = -1 [default] colvars: The force constant for colvar "d" will be rescaled to 250 according to the specified width (0.2). colvars: ---------------------------------------------------------------------- @@ -114,18 +118,38 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 2 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: +colvars: SUMMARY OF COLVARS FEATURES USED SO FAR AND THEIR CITATIONS: +colvars: +colvars: - Colvars module: +colvars: - Colvars-NAMD interface: +colvars: - Harmonic colvar bias implementation: +colvars: - Optimal rotation via flexible fitting: +colvars: - distance colvar component: +colvars: Fiorin2013 https://doi.org/10.1080/00268976.2013.813594 +colvars: +colvars: - ABF colvar bias implementation: +colvars: - Internal-forces free energy estimator: +colvars: Henin2010 https://doi.org/10.1021/ct9004432 +colvars: +colvars: - NAMD engine: +colvars: - Scalable center-of-mass computation (NAMD): +colvars: Phillips2020 https://doi.org/10.1063/5.0014475 +colvars: +colvars: updating target temperature (T = 300 K). +colvars: Current simulation parameters: initial step = 0, integration timestep = 0.5 +colvars: Updating atomic parameters (masses, charges, etc). colvars: Re-initialized atom group for variable "d":0/0. 1 atoms: total mass = 12.011, total charge = -0.1. colvars: Re-initialized atom group for variable "d":0/1. 1 atoms: total mass = 12.011, total charge = 0.07. colvars: The final output state file will be "test.colvars.state". -colvars: Opening trajectory file "test.colvars.traj". -colvars: Redefining the Tcl "cv" command to the new script interface. -colvars: Prepared sample and gradient buffers at step 0. colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 2 atoms in total). -colvars: Re-initialized atom group for variable "d":0/0. 1 atoms: total mass = 12.011, total charge = -0.1. -colvars: Re-initialized atom group for variable "d":0/1. 1 atoms: total mass = 12.011, total charge = 0.07. -colvars: The final output state file will be "test.colvars.state". +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 300 K). +colvars: Current simulation parameters: initial step = 0, integration timestep = 0.5 +colvars: Updating atomic parameters (masses, charges, etc). colvars: Saving collective variables state to "test.colvars.state". colvars: ---------------------------------------------------------------------- colvars: Reading new configuration: @@ -135,6 +159,7 @@ colvars: # colvarsTrajFrequency = 2 [default] colvars: # colvarsRestartFrequency = 0 [default] colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Collective variables initialized, 1 in total. colvars: ---------------------------------------------------------------------- @@ -151,6 +176,7 @@ colvars: # centers = { 13 } colvars: # targetCenters = { 13 } [default] colvars: # outputCenters = off [default] colvars: # forceConstant = 1 +colvars: # decoupling = off [default] colvars: # targetForceConstant = -1 [default] colvars: The force constant for colvar "d" will be rescaled to 25 according to the specified width (0.2). colvars: ---------------------------------------------------------------------- @@ -161,273 +187,296 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 2 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 300 K). colvars: Saving collective variables state to "test.colvars.state". -colvars: Features of "colvar d" ON/OFF (refcount) -colvars: - active ON (4) -colvars: - awake ON (0) -colvars: - gradient ON (3) -colvars: - collect_gradient OFF (0) -colvars: - velocity_from_finite_differences OFF (0) -colvars: - total_force ON (3) -colvars: - total_force_calculation ON (1) -colvars: - subtract_applied_force_from_total_force ON (0) -colvars: - Jacobian_derivative ON (1) -colvars: - hide_Jacobian_force OFF (0) -colvars: - extended_Lagrangian OFF (0) -colvars: - Langevin_dynamics OFF (0) -colvars: - output_energy OFF (0) -colvars: - output_value ON (0) -colvars: - output_velocity OFF (0) -colvars: - output_applied_force ON (0) -colvars: - output_total_force ON (0) -colvars: - lower_boundary ON (1) -colvars: - upper_boundary ON (1) -colvars: - hard_lower_boundary OFF (0) -colvars: - hard_upper_boundary OFF (0) -colvars: - reflecting_lower_boundary OFF (0) -colvars: - reflecting_upper_boundary OFF (0) -colvars: - grid ON (0) -colvars: - running_average OFF (0) -colvars: - correlation_function OFF (0) -colvars: - scripted OFF (0) -colvars: - custom_function OFF (0) -colvars: - periodic OFF (0) -colvars: - single_component ON (0) -colvars: - scalar ON (6) -colvars: - linear ON (2) -colvars: - homogeneous ON (0) -colvars: - multiple_timestep OFF (0) +colvars: Features of "colvar d" (refcount) +colvars: - ON active (3) +colvars: - awake +colvars: - ON gradient (3) +colvars: - collect_gradient +colvars: - collect_atom_ids +colvars: - velocity_from_finite_differences +colvars: - ON total_force (3) +colvars: - ON total_force_calculation (1) +colvars: - ON subtract_applied_force_from_total_force +colvars: - ON Jacobian_derivative (1) +colvars: - hide_Jacobian_force +colvars: - extended_Lagrangian +colvars: - external +colvars: - Langevin_dynamics +colvars: - output_energy +colvars: - ON output_value +colvars: - output_velocity +colvars: - ON output_applied_force +colvars: - ON output_total_force +colvars: - ON lower_boundary +colvars: - ON upper_boundary +colvars: - hard_lower_boundary +colvars: - hard_upper_boundary +colvars: - reflecting_lower_boundary +colvars: - reflecting_upper_boundary +colvars: - ON grid +colvars: - running_average +colvars: - correlation_function +colvars: - scripted +colvars: - custom_function +colvars: - periodic +colvars: - ON single_component +colvars: - ON scalar (7) +colvars: - ON linear (2) +colvars: - ON homogeneous +colvars: - multiple_timestep colvars: * child 1 -colvars: Features of "cvc distance0001" ON/OFF (refcount) -colvars: - active ON (0) -colvars: - scalar OFF (0) -colvars: - periodic OFF (0) -colvars: - defined_width OFF (0) -colvars: - defined_lower_boundary ON (0) -colvars: - defined_upper_boundary OFF (0) -colvars: - gradient ON (2) -colvars: - explicit_gradient ON (0) -colvars: - inverse_gradient ON (2) -colvars: - Jacobian_derivative ON (1) -colvars: - total_force_from_one_group OFF (0) -colvars: - debug_gradient OFF (0) -colvars: - use_minimum-image_with_PBCs ON (0) -colvars: - function_of_centers_of_mass ON (1) -colvars: - scalable_calculation ON (0) -colvars: - scalable_calculation_of_centers_of_mass ON (1) +colvars: Features of "cvc distance0001 of type "distance"" (refcount) +colvars: - ON active +colvars: - scalar +colvars: - periodic +colvars: - defined_width +colvars: - ON defined_lower_boundary +colvars: - defined_upper_boundary +colvars: - ON gradient (2) +colvars: - explicit_gradient +colvars: - ON inverse_gradient (2) +colvars: - ON Jacobian_derivative (1) +colvars: - total_force_from_one_group +colvars: - debug_gradient +colvars: - ON use_minimum-image_with_PBCs +colvars: - ON function_of_centers_of_mass (1) +colvars: - ON scalable_calculation (1) +colvars: - ON scalable_calculation_of_centers_of_mass +colvars: - collect_atom_ids colvars: * child 1 -colvars: Features of "atom group group1" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) +colvars: Features of "atom group group1" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids colvars: * child 2 -colvars: Features of "atom group group2" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) -colvars: Features of "bias harmonic1" ON/OFF (refcount) -colvars: - active ON (1) -colvars: - awake ON (0) -colvars: - step_zero_data OFF (0) -colvars: - apply_force ON (0) -colvars: - bypass_extended_Lagrangian_coordinates OFF (0) -colvars: - obtain_total_force OFF (0) -colvars: - output_accumulated_work OFF (0) -colvars: - history_dependent OFF (0) -colvars: - time_dependent OFF (0) -colvars: - require_scalar_variables OFF (0) -colvars: - calculate_a_PMF OFF (0) -colvars: - calculate_TI_samples OFF (0) -colvars: - write_TI_samples_ OFF (0) -colvars: - write_TI_PMF OFF (0) +colvars: Features of "atom group group2" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids +colvars: Features of "bias harmonic1" (refcount) +colvars: - ON active +colvars: - awake +colvars: - step_zero_data +colvars: - ON apply_force +colvars: - bypass_extended_Lagrangian_coordinates +colvars: - obtain_total_force +colvars: - output_accumulated_work +colvars: - history_dependent +colvars: - time_dependent +colvars: - require_scalar_variables +colvars: - calculate_a_PMF +colvars: - calculate_TI_samples +colvars: - write_TI_samples_ +colvars: - write_TI_PMF +colvars: - scale_biasing_force colvars: * child 1 -colvars: Features of "colvar d" ON/OFF (refcount) -colvars: - active ON (4) -colvars: - awake ON (0) -colvars: - gradient ON (3) -colvars: - collect_gradient OFF (0) -colvars: - velocity_from_finite_differences OFF (0) -colvars: - total_force ON (3) -colvars: - total_force_calculation ON (1) -colvars: - subtract_applied_force_from_total_force ON (0) -colvars: - Jacobian_derivative ON (1) -colvars: - hide_Jacobian_force OFF (0) -colvars: - extended_Lagrangian OFF (0) -colvars: - Langevin_dynamics OFF (0) -colvars: - output_energy OFF (0) -colvars: - output_value ON (0) -colvars: - output_velocity OFF (0) -colvars: - output_applied_force ON (0) -colvars: - output_total_force ON (0) -colvars: - lower_boundary ON (1) -colvars: - upper_boundary ON (1) -colvars: - hard_lower_boundary OFF (0) -colvars: - hard_upper_boundary OFF (0) -colvars: - reflecting_lower_boundary OFF (0) -colvars: - reflecting_upper_boundary OFF (0) -colvars: - grid ON (0) -colvars: - running_average OFF (0) -colvars: - correlation_function OFF (0) -colvars: - scripted OFF (0) -colvars: - custom_function OFF (0) -colvars: - periodic OFF (0) -colvars: - single_component ON (0) -colvars: - scalar ON (6) -colvars: - linear ON (2) -colvars: - homogeneous ON (0) -colvars: - multiple_timestep OFF (0) +colvars: Features of "colvar d" (refcount) +colvars: - ON active (3) +colvars: - awake +colvars: - ON gradient (3) +colvars: - collect_gradient +colvars: - collect_atom_ids +colvars: - velocity_from_finite_differences +colvars: - ON total_force (3) +colvars: - ON total_force_calculation (1) +colvars: - ON subtract_applied_force_from_total_force +colvars: - ON Jacobian_derivative (1) +colvars: - hide_Jacobian_force +colvars: - extended_Lagrangian +colvars: - external +colvars: - Langevin_dynamics +colvars: - output_energy +colvars: - ON output_value +colvars: - output_velocity +colvars: - ON output_applied_force +colvars: - ON output_total_force +colvars: - ON lower_boundary +colvars: - ON upper_boundary +colvars: - hard_lower_boundary +colvars: - hard_upper_boundary +colvars: - reflecting_lower_boundary +colvars: - reflecting_upper_boundary +colvars: - ON grid +colvars: - running_average +colvars: - correlation_function +colvars: - scripted +colvars: - custom_function +colvars: - periodic +colvars: - ON single_component +colvars: - ON scalar (7) +colvars: - ON linear (2) +colvars: - ON homogeneous +colvars: - multiple_timestep colvars: * child 1 -colvars: Features of "cvc distance0001" ON/OFF (refcount) -colvars: - active ON (0) -colvars: - scalar OFF (0) -colvars: - periodic OFF (0) -colvars: - defined_width OFF (0) -colvars: - defined_lower_boundary ON (0) -colvars: - defined_upper_boundary OFF (0) -colvars: - gradient ON (2) -colvars: - explicit_gradient ON (0) -colvars: - inverse_gradient ON (2) -colvars: - Jacobian_derivative ON (1) -colvars: - total_force_from_one_group OFF (0) -colvars: - debug_gradient OFF (0) -colvars: - use_minimum-image_with_PBCs ON (0) -colvars: - function_of_centers_of_mass ON (1) -colvars: - scalable_calculation ON (0) -colvars: - scalable_calculation_of_centers_of_mass ON (1) +colvars: Features of "cvc distance0001 of type "distance"" (refcount) +colvars: - ON active +colvars: - scalar +colvars: - periodic +colvars: - defined_width +colvars: - ON defined_lower_boundary +colvars: - defined_upper_boundary +colvars: - ON gradient (2) +colvars: - explicit_gradient +colvars: - ON inverse_gradient (2) +colvars: - ON Jacobian_derivative (1) +colvars: - total_force_from_one_group +colvars: - debug_gradient +colvars: - ON use_minimum-image_with_PBCs +colvars: - ON function_of_centers_of_mass (1) +colvars: - ON scalable_calculation (1) +colvars: - ON scalable_calculation_of_centers_of_mass +colvars: - collect_atom_ids colvars: * child 1 -colvars: Features of "atom group group1" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) +colvars: Features of "atom group group1" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids colvars: * child 2 -colvars: Features of "atom group group2" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) +colvars: Features of "atom group group2" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids colvars: Saving collective variables state to "test.colvars.state". -colvars: Features of "bias harm2" ON/OFF (refcount) -colvars: - active ON (1) -colvars: - awake ON (0) -colvars: - step_zero_data OFF (0) -colvars: - apply_force ON (0) -colvars: - bypass_extended_Lagrangian_coordinates OFF (0) -colvars: - obtain_total_force OFF (0) -colvars: - output_accumulated_work OFF (0) -colvars: - history_dependent OFF (0) -colvars: - time_dependent OFF (0) -colvars: - require_scalar_variables OFF (0) -colvars: - calculate_a_PMF OFF (0) -colvars: - calculate_TI_samples OFF (0) -colvars: - write_TI_samples_ OFF (0) -colvars: - write_TI_PMF OFF (0) +colvars: Features of "bias harm2" (refcount) +colvars: - ON active +colvars: - awake +colvars: - step_zero_data +colvars: - ON apply_force +colvars: - bypass_extended_Lagrangian_coordinates +colvars: - obtain_total_force +colvars: - output_accumulated_work +colvars: - history_dependent +colvars: - time_dependent +colvars: - require_scalar_variables +colvars: - calculate_a_PMF +colvars: - calculate_TI_samples +colvars: - write_TI_samples_ +colvars: - write_TI_PMF +colvars: - scale_biasing_force colvars: * child 1 -colvars: Features of "colvar d" ON/OFF (refcount) -colvars: - active ON (3) -colvars: - awake ON (0) -colvars: - gradient ON (2) -colvars: - collect_gradient OFF (0) -colvars: - velocity_from_finite_differences OFF (0) -colvars: - total_force ON (3) -colvars: - total_force_calculation ON (1) -colvars: - subtract_applied_force_from_total_force ON (0) -colvars: - Jacobian_derivative ON (1) -colvars: - hide_Jacobian_force OFF (0) -colvars: - extended_Lagrangian OFF (0) -colvars: - Langevin_dynamics OFF (0) -colvars: - output_energy OFF (0) -colvars: - output_value ON (0) -colvars: - output_velocity OFF (0) -colvars: - output_applied_force ON (0) -colvars: - output_total_force ON (0) -colvars: - lower_boundary ON (1) -colvars: - upper_boundary ON (1) -colvars: - hard_lower_boundary OFF (0) -colvars: - hard_upper_boundary OFF (0) -colvars: - reflecting_lower_boundary OFF (0) -colvars: - reflecting_upper_boundary OFF (0) -colvars: - grid ON (0) -colvars: - running_average OFF (0) -colvars: - correlation_function OFF (0) -colvars: - scripted OFF (0) -colvars: - custom_function OFF (0) -colvars: - periodic OFF (0) -colvars: - single_component ON (0) -colvars: - scalar ON (6) -colvars: - linear ON (2) -colvars: - homogeneous ON (0) -colvars: - multiple_timestep OFF (0) +colvars: Features of "colvar d" (refcount) +colvars: - ON active (2) +colvars: - awake +colvars: - ON gradient (2) +colvars: - collect_gradient +colvars: - collect_atom_ids +colvars: - velocity_from_finite_differences +colvars: - ON total_force (3) +colvars: - ON total_force_calculation (1) +colvars: - ON subtract_applied_force_from_total_force +colvars: - ON Jacobian_derivative (1) +colvars: - hide_Jacobian_force +colvars: - extended_Lagrangian +colvars: - external +colvars: - Langevin_dynamics +colvars: - output_energy +colvars: - ON output_value +colvars: - output_velocity +colvars: - ON output_applied_force +colvars: - ON output_total_force +colvars: - ON lower_boundary +colvars: - ON upper_boundary +colvars: - hard_lower_boundary +colvars: - hard_upper_boundary +colvars: - reflecting_lower_boundary +colvars: - reflecting_upper_boundary +colvars: - ON grid +colvars: - running_average +colvars: - correlation_function +colvars: - scripted +colvars: - custom_function +colvars: - periodic +colvars: - ON single_component +colvars: - ON scalar (7) +colvars: - ON linear (2) +colvars: - ON homogeneous +colvars: - multiple_timestep colvars: * child 1 -colvars: Features of "cvc distance0001" ON/OFF (refcount) -colvars: - active ON (0) -colvars: - scalar OFF (0) -colvars: - periodic OFF (0) -colvars: - defined_width OFF (0) -colvars: - defined_lower_boundary ON (0) -colvars: - defined_upper_boundary OFF (0) -colvars: - gradient ON (2) -colvars: - explicit_gradient ON (0) -colvars: - inverse_gradient ON (2) -colvars: - Jacobian_derivative ON (1) -colvars: - total_force_from_one_group OFF (0) -colvars: - debug_gradient OFF (0) -colvars: - use_minimum-image_with_PBCs ON (0) -colvars: - function_of_centers_of_mass ON (1) -colvars: - scalable_calculation ON (0) -colvars: - scalable_calculation_of_centers_of_mass ON (1) +colvars: Features of "cvc distance0001 of type "distance"" (refcount) +colvars: - ON active +colvars: - scalar +colvars: - periodic +colvars: - defined_width +colvars: - ON defined_lower_boundary +colvars: - defined_upper_boundary +colvars: - ON gradient (2) +colvars: - explicit_gradient +colvars: - ON inverse_gradient (2) +colvars: - ON Jacobian_derivative (1) +colvars: - total_force_from_one_group +colvars: - debug_gradient +colvars: - ON use_minimum-image_with_PBCs +colvars: - ON function_of_centers_of_mass (1) +colvars: - ON scalable_calculation (1) +colvars: - ON scalable_calculation_of_centers_of_mass +colvars: - collect_atom_ids colvars: * child 1 -colvars: Features of "atom group group1" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) +colvars: Features of "atom group group1" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids colvars: * child 2 -colvars: Features of "atom group group2" ON/OFF (refcount) -colvars: - active OFF (0) -colvars: - translational_fit OFF (0) -colvars: - rotational_fit OFF (0) -colvars: - fitting_group OFF (0) -colvars: - explicit_atom_gradient ON (1) -colvars: - fit_gradients OFF (0) -colvars: - atomic_forces OFF (0) -colvars: - scalable_group ON (0) -colvars: - scalable_group_center_of_mass ON (1) -colvars: The final output state file will be "test.colvars.state". +colvars: Features of "atom group group2" (refcount) +colvars: - active +colvars: - center_to_reference +colvars: - center_to_origin +colvars: - rotate_to_origin +colvars: - fitting_group +colvars: - explicit_atom_gradient +colvars: - fit_gradients +colvars: - atomic_forces +colvars: - ON scalable_group (1) +colvars: - ON scalable_group_center_of_mass +colvars: - collect_atom_ids colvars: Saving collective variables state to "test.colvars.state". -colvars: The final output state file will be "test.colvars.state". colvars: Saving collective variables state to "test.colvars.state". colvars: ---------------------------------------------------------------------- -colvars: Restarting from file "test.colvars.state". +colvars: Loading state from text file "test.colvars.state". colvars: Restarting collective variable "d" from value: 11.9615 -colvars: Restarting abf bias "abf1" from step number 60. +colvars: Restarted abf bias "abf1" with step number 60. colvars: ---------------------------------------------------------------------- colvars: Resetting the Collective Variables module. colvars: ---------------------------------------------------------------------- @@ -438,6 +487,7 @@ colvars: # colvarsTrajFrequency = 2 colvars: # colvarsRestartFrequency = 0 [default] colvars: # scriptedColvarForces = off [default] colvars: # scriptingAfterBiases = off [default] +colvars: # sourceTclFile = "" [default] colvars: ---------------------------------------------------------------------- colvars: Initializing a new collective variable. colvars: # name = "d" @@ -451,32 +501,32 @@ colvars: # forceNoPBC = off [default] colvars: # scalable = on [default] colvars: Initializing atom group "group1". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "" [default] colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "group1". colvars: # printAtomIDs = off [default] -colvars: Atom group "group1" defined with 1 atoms requested: total mass = 12.011, total charge = -0.1. +colvars: Atom group "group1" defined with 1 atoms requested. colvars: Initializing atom group "group2". colvars: # name = "" [default] -colvars: # centerReference = off [default] -colvars: # rotateReference = off [default] +colvars: # centerToOrigin = off [default] +colvars: # centerToReference = off [default] +colvars: # rotateToReference = off [default] colvars: # atomsOfGroup = "" [default] colvars: # indexGroup = "" [default] colvars: # psfSegID = [default] colvars: # atomsFile = "" [default] colvars: # dummyAtom = ( 0 , 0 , 0 ) [default] -colvars: # enableForces = on [default] colvars: # enableFitGradients = on [default] colvars: Enabling scalable calculation for group "group2". colvars: # printAtomIDs = off [default] -colvars: Atom group "group2" defined with 1 atoms requested: total mass = 12.011, total charge = 0.07. +colvars: Atom group "group2" defined with 1 atoms requested. colvars: # oneSiteSystemForce = off [default] colvars: # oneSiteTotalForce = off [default] colvars: All components initialized. @@ -506,13 +556,14 @@ colvars: # outputEnergy = off [default] colvars: # outputFreq = 60 colvars: # timeStepFactor = 1 [default] colvars: # applyBias = on [default] -colvars: # updateBias = on [default] colvars: # hideJacobian = off [default] colvars: Jacobian (geometric) forces will be included in reported free energy gradients. colvars: # fullSamples = 0 +colvars: # minSamples = 0 [default] colvars: # inputPrefix = [default] colvars: # historyFreq = 0 [default] colvars: # shared = off [default] +colvars: # updateBias = on [default] colvars: # maxForce = [default] colvars: # integrate = on [default] colvars: Finished ABF setup. @@ -530,6 +581,7 @@ colvars: # centers = { 12 } colvars: # targetCenters = { 12 } [default] colvars: # outputCenters = off [default] colvars: # forceConstant = 10 +colvars: # decoupling = off [default] colvars: # targetForceConstant = -1 [default] colvars: The force constant for colvar "d" will be rescaled to 250 according to the specified width (0.2). colvars: ---------------------------------------------------------------------- @@ -540,5 +592,7 @@ colvars: ---------------------------------------------------------------------- colvars: Updating NAMD interface: colvars: updating atomic data (0 atoms). colvars: updating group data (2 scalable groups, 2 atoms in total). +colvars: updating grid object data (0 grid objects in total). +colvars: updating target temperature (T = 300 K). colvars: Warning: before deleting colvar d, deleting related bias harmonic1 colvars: Warning: before deleting colvar d, deleting related bias abf1 diff --git a/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.traj b/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.traj index f029f1ae9..d8c47efc5 100644 --- a/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.traj +++ b/namd/tests/library/005_Tcl_script/AutoDiff/test.colvars.traj @@ -1,38 +1,38 @@ # step d ft_d fa_d E_harmonic1 - 0 0.00000000000000e+00 0.00000000000000e+00 3.00000000000000e+03 1.80000000000000e+04 + 0 0.00000000000000e+00 0.00000000000000e+00 3.00000000000000e+03 1.80000000000000e+04 # step d ft_d fa_d E_harmonic1 - 0 1.42467449615693e+01 0.00000000000000e+00 -5.61686240392322e+02 6.30982865292123e+02 - 2 1.42282559728026e+01 5.32933601891204e+00 -5.59529927203773e+02 6.20640585041317e+02 - 4 1.41727662460744e+01 3.36254334374946e+01 -5.43191561518601e+02 5.90114145010033e+02 - 6 1.40831323320990e+01 8.13697354080271e+01 -5.89087882343440e+02 5.42430039129534e+02 - 8 1.39609835871970e+01 1.43683712693643e+02 -4.90245896799256e+02 4.80682078657013e+02 - 10 1.38136915830229e+01 2.11747525548495e+02 -6.48260269460101e+02 4.11184644791030e+02 - 12 1.36437184346326e+01 2.81640059905252e+02 -6.92569668563411e+02 3.37726286543895e+02 - 14 1.34588898265340e+01 3.45729335512429e+02 -7.10451792145928e+02 2.66044940745548e+02 - 16 1.32645132385777e+01 4.00985013276056e+02 -7.17113322920477e+02 1.99874216317278e+02 - 18 1.30658878095430e+01 4.46271310459046e+02 -7.12743262844788e+02 1.42014602816538e+02 - 20 1.28684277602171e+01 4.80580457979151e+02 -6.97687398033433e+02 9.42708468394674e+01 + 0 1.42467449615693e+01 0.00000000000000e+00 -5.61686240392322e+02 6.30982865292123e+02 + 2 1.42282559728026e+01 5.32933601659170e+00 -5.59529927201590e+02 6.20640585041297e+02 + 4 1.41727662460742e+01 3.36254334328178e+01 -5.43191561518562e+02 5.90114145009947e+02 + 6 1.40831323320984e+01 8.13697353990826e+01 -5.89087882335502e+02 5.42430039129236e+02 + 8 1.39609835871958e+01 1.43683712679872e+02 -4.90245896798950e+02 4.80682078656413e+02 + 10 1.38136915830204e+01 2.11747525531865e+02 -6.48260269443340e+02 4.11184644789902e+02 + 12 1.36437184346286e+01 2.81640059890150e+02 -6.92569668547312e+02 3.37726286542257e+02 + 14 1.34588898265283e+01 3.45729335505761e+02 -7.10451792137843e+02 2.66044940743482e+02 + 16 1.32645132385707e+01 4.00985013285219e+02 -7.17113322927902e+02 1.99874216315080e+02 + 18 1.30658878095359e+01 4.46271310489825e+02 -7.12743262873798e+02 1.42014602814653e+02 + 20 1.28684277602119e+01 4.80580458032557e+02 -6.97687398085530e+02 9.42708468383301e+01 # step d ft_d fa_d E_harmonic1 E_harm2 - 20 1.28684277602171e+01 4.91888943151427e+02 -7.00052334624999e+02 9.42708468394674e+01 2.16390678518517e-01 - 22 1.26751474320475e+01 5.03084201909373e+02 -6.63749745722433e+02 5.69780068750397e+01 1.31911488631676e+00 - 24 1.24888492973110e+01 5.15335038327460e+02 -6.24768595087998e+02 2.98717044351878e+01 3.26593801074265e+00 - 26 1.23139243019631e+01 5.23015189033541e+02 -5.84344372073400e+02 1.23185584203791e+01 5.88374829295988e+00 - 28 1.21466432033769e+01 5.30527139707869e+02 -1.53268809286588e+01 2.68802863708136e+00 9.10272277928461e+00 - 30 1.19993450537562e+01 5.38695710361251e+02 2.51801102170466e+01 5.36193227844337e-05 1.25163790180274e+01 - 32 1.18850444890963e+01 5.34994463383540e+02 5.66127654985307e+01 1.65184618589271e+00 1.55390723911830e+01 - 34 1.18115458381165e+01 5.01682267555293e+02 7.68248945179519e+01 4.43937139139973e+00 1.76552911862265e+01 - 36 1.17769967622096e+01 4.41508793313709e+02 8.63258903923681e+01 6.21630550812685e+00 1.86967114955734e+01 - 38 1.17769542139854e+01 3.61640408006584e+02 8.63375911540067e+01 6.21867783235716e+00 1.86980124336000e+01 - 40 1.18061287066284e+01 2.69866689808628e+02 7.83146056771794e+01 4.69825979919525e+00 1.78166083142086e+01 + 20 1.28684277602119e+01 4.91888943214004e+02 -6.94398092090827e+02 9.42708468383301e+01 2.16390678535748e-01 + 22 1.26752449413069e+01 5.03055897841674e+02 -6.63748256701074e+02 5.69944663450714e+01 1.31832310183441e+00 + 24 1.24890388622314e+01 5.15261138484907e+02 -6.24746825598548e+02 2.98948760965754e+01 3.26351605387201e+00 + 26 1.23141975756229e+01 5.22905582834307e+02 -5.84309916130603e+02 1.23400145659133e+01 5.87906206601890e+00 + 28 1.21469912236785e+01 5.30384209752227e+02 -1.54225865115908e+01 2.70080247981329e+00 9.09529965601853e+00 + 30 1.19997542475315e+01 5.38523525126186e+02 2.50675819288300e+01 7.54928447006639e-06 1.25061445666403e+01 + 32 1.18854965342920e+01 5.34804689652633e+02 5.64884530697052e+01 1.63888045739343e+00 1.55264746884398e+01 + 34 1.18120184515516e+01 5.01493410173479e+02 7.66949258233203e+01 4.41713281963403e+00 1.76412519931743e+01 + 36 1.17774680364548e+01 4.41334033478252e+02 8.61962899749342e+01 6.19005934991112e+00 1.86823050236215e+01 + 38 1.17774037378490e+01 3.61487108034449e+02 8.62139720915345e+01 6.19363699045151e+00 1.86842702528210e+01 + 40 1.18065383593429e+01 2.69739676343623e+02 7.82019511807010e+01 4.67842580071685e+00 1.78043835964990e+01 # step d ft_d fa_d E_harm2 - 40 1.18061287066284e+01 2.21949532217238e+02 2.98467823342890e+01 1.78166083142086e+01 - 42 1.18555777458287e+01 1.74464230745251e+02 2.86105563542814e+01 1.63712786980303e+01 - 44 1.19176677646512e+01 8.22937404700844e+01 2.70583058837192e+01 1.46430383459383e+01 - 46 1.19851813910175e+01 -6.16868563783110e+00 2.53704652245623e+01 1.28732101142145e+01 - 48 1.20415309208210e+01 -8.97664352869559e+01 -3.04621864055419e+02 1.14832871967769e+01 - 50 1.20655400380907e+01 -1.50412424518311e+02 -1.19070572741895e+02 1.09151927551443e+01 - 52 1.20686235841253e+01 -1.88148861003765e+02 -2.70255157752960e+01 1.08432753505955e+01 - 54 1.20571466466005e+01 -2.05067811286631e+02 2.93687990112175e+01 1.11121555752091e+01 - 56 1.20332538278214e+01 -2.00399869950163e+02 6.58114558084562e+01 1.16824770177739e+01 - 58 1.20000609172601e+01 -1.75170570145191e+02 8.83273231404457e+01 1.24984771148846e+01 - 60 1.19614725548519e+01 -1.36757725550194e+02 2.59631861287023e+01 1.34817406790728e+01 + 40 1.18065383593429e+01 2.21838011341288e+02 2.98365410164274e+01 1.78043835964990e+01 + 42 1.18559358360142e+01 1.74370020065540e+02 2.86016040996450e+01 1.63610351414566e+01 + 44 1.19179690389992e+01 8.22420089611619e+01 2.70507740250190e+01 1.46348875070528e+01 + 46 1.19854241580116e+01 -6.16498136864664e+00 2.53643960497090e+01 1.28670517393298e+01 + 48 1.20417205707688e+01 -8.97018803497429e+01 -3.04537649336825e+02 1.14787433060950e+01 + 50 1.20656896196389e+01 -1.50302554047029e+02 -1.19060978418300e+02 1.09116985856314e+01 + 52 1.20687425752562e+01 -1.88023807239229e+02 -2.70542194599442e+01 1.08405048892546e+01 + 54 1.20572415396545e+01 -2.04963680940880e+02 2.93216002704265e+01 1.11099189319118e+01 + 56 1.20333270708426e+01 -2.00351628760287e+02 6.57613077254422e+01 1.16807068995712e+01 + 58 1.20001102284550e+01 -1.75192284978711e+02 8.82857871369106e+01 1.24972444405037e+01 + 60 1.19614929506509e+01 -1.36804624327918e+02 2.59626762337285e+01 1.34812111443482e+01 diff --git a/namd/tests/library/005_Tcl_script/AutoDiff/test.grad b/namd/tests/library/005_Tcl_script/AutoDiff/test.grad index 2a30d2884..9d0cba4b1 100644 --- a/namd/tests/library/005_Tcl_script/AutoDiff/test.grad +++ b/namd/tests/library/005_Tcl_script/AutoDiff/test.grad @@ -1,18 +1,18 @@ # 1 # 1.20000000000000e+01 2.00000000000000e-01 100 0 - 1.21000000000000e+01 6.99836008077528e+01 - 1.23000000000000e+01 -5.26375313396654e+02 - 1.25000000000000e+01 -5.16708951944646e+02 - 1.27000000000000e+01 -5.05977066222890e+02 - 1.29000000000000e+01 -4.88119448094002e+02 - 1.31000000000000e+01 -4.54995001210363e+02 - 1.33000000000000e+01 -4.12314411188385e+02 - 1.35000000000000e+01 -3.59715820941386e+02 - 1.37000000000000e+01 -2.97680997044662e+02 - 1.39000000000000e+01 -2.12126120236132e+02 - 1.41000000000000e+01 -9.78647225049371e+01 - 1.43000000000000e+01 -1.38644228524206e+01 + 1.21000000000000e+01 6.99496988696531e+01 + 1.23000000000000e+01 -5.26249052199959e+02 + 1.25000000000000e+01 -5.16625922714804e+02 + 1.27000000000000e+01 -5.05936450367739e+02 + 1.29000000000000e+01 -4.86234700623280e+02 + 1.31000000000000e+01 -4.54995001246909e+02 + 1.33000000000000e+01 -4.12314411202659e+02 + 1.35000000000000e+01 -3.59715820938218e+02 + 1.37000000000000e+01 -2.97680997031215e+02 + 1.39000000000000e+01 -2.12126120219849e+02 + 1.41000000000000e+01 -9.78647224947473e+01 + 1.43000000000000e+01 -1.38644228493631e+01 1.45000000000000e+01 0.00000000000000e+00 1.47000000000000e+01 0.00000000000000e+00 1.49000000000000e+01 0.00000000000000e+00 diff --git a/namd/tests/library/005_Tcl_script/AutoDiff/test.pmf b/namd/tests/library/005_Tcl_script/AutoDiff/test.pmf index edf8b154a..7d017e393 100644 --- a/namd/tests/library/005_Tcl_script/AutoDiff/test.pmf +++ b/namd/tests/library/005_Tcl_script/AutoDiff/test.pmf @@ -1,18 +1,18 @@ # 1 # 1.19000000000000e+01 2.00000000000000e-01 101 0 - 1.20000000000000e+01 7.63151734965745e+02 - 1.22000000000000e+01 7.77148455127296e+02 - 1.24000000000000e+01 6.71873392447965e+02 - 1.26000000000000e+01 5.68531602059036e+02 - 1.28000000000000e+01 4.67336188814458e+02 - 1.30000000000000e+01 3.69712299195657e+02 - 1.32000000000000e+01 2.78713298953585e+02 - 1.34000000000000e+01 1.96250416715908e+02 - 1.36000000000000e+01 1.24307252527630e+02 - 1.38000000000000e+01 6.47710531186980e+01 - 1.40000000000000e+01 2.23458290714716e+01 - 1.42000000000000e+01 2.77288457048417e+00 + 1.20000000000000e+01 7.62731584603818e+02 + 1.22000000000000e+01 7.76721524377749e+02 + 1.24000000000000e+01 6.71471713937757e+02 + 1.26000000000000e+01 5.68146529394796e+02 + 1.28000000000000e+01 4.66959239321248e+02 + 1.30000000000000e+01 3.69712299196592e+02 + 1.32000000000000e+01 2.78713298947211e+02 + 1.34000000000000e+01 1.96250416706679e+02 + 1.36000000000000e+01 1.24307252519035e+02 + 1.38000000000000e+01 6.47710531127921e+01 + 1.40000000000000e+01 2.23458290688221e+01 + 1.42000000000000e+01 2.77288456987264e+00 1.44000000000000e+01 0.00000000000000e+00 1.46000000000000e+01 0.00000000000000e+00 1.48000000000000e+01 0.00000000000000e+00 diff --git a/namd/tests/library/multiple_walker_abf/abf_12-15_less.colvars b/namd/tests/library/multiple_walker_abf/abf_12-15_less.colvars index 84d1a9488..6b7e392a7 100644 --- a/namd/tests/library/multiple_walker_abf/abf_12-15_less.colvars +++ b/namd/tests/library/multiple_walker_abf/abf_12-15_less.colvars @@ -28,7 +28,6 @@ atomNumbers { 10 } } abf { - name abf colvars endToEnd fullSamples 500 historyFreq 20000 diff --git a/namd/tests/library/multiple_walker_abf/mwABF_selection.tcl b/namd/tests/library/multiple_walker_abf/mwABF_selection.tcl new file mode 120000 index 000000000..3c728e9fa --- /dev/null +++ b/namd/tests/library/multiple_walker_abf/mwABF_selection.tcl @@ -0,0 +1 @@ +../../../../colvartools/mwABF_selection.tcl \ No newline at end of file diff --git a/namd/tests/library/multiple_walker_abf/mw_independent.namd b/namd/tests/library/multiple_walker_abf/mw_independent.namd index 0b14a5689..05f24644e 100644 --- a/namd/tests/library/multiple_walker_abf/mw_independent.namd +++ b/namd/tests/library/multiple_walker_abf/mw_independent.namd @@ -1,7 +1,7 @@ # set: thisName lastName run temp reinitTemp outputName out_mw_independent.[myReplica] print "REPLICA [myReplica]" -seed [expr {[myReplica]*[myReplica]*[myReplica]*2803 + 99377}] +seed [expr {[myReplica]*[myReplica]*[myReplica]*2803 + 1}] set temp 300 numSteps 40000 diff --git a/namd/tests/library/multiple_walker_abf/mw_selection.namd b/namd/tests/library/multiple_walker_abf/mw_selection.namd index bcf86322a..e045be9f8 100644 --- a/namd/tests/library/multiple_walker_abf/mw_selection.namd +++ b/namd/tests/library/multiple_walker_abf/mw_selection.namd @@ -1,11 +1,11 @@ # set: thisName lastName run temp reinitTemp -outputName out_mw_selection.[myReplica] +outputName out_mw_selection.[myReplica] print "REPLICA [myReplica]" seed [expr {[myReplica]*[myReplica]*[myReplica]*20233 + 20219}] set temp 300 -structure ../Common/da.psf -coordinates ../Common/da.min.pdb +structure ../Common/da.psf +coordinates ../Common/da.min.pdb temperature $temp cellBasisVector1 200 0 0 cellBasisVector2 0 200 0 @@ -15,7 +15,7 @@ cellBasisVector3 0 0 200 langevin on langevinTemp $temp langevinHydrogen off -langevinDamping 1 +langevinDamping 1 # parameters parameters ../Common/par_all22_prot.inp @@ -38,157 +38,29 @@ stepsPerCycle 10 # output binaryOutput yes binaryRestart yes -wrapAll yes -wrapNearest yes +wrapAll yes +wrapNearest yes comMotion yes outputEnergies 1000 -#outputPressure 1000 +#outputPressure 1000 outputTiming 1000 #xstFreq 1000 dcdFreq 5000 restartFreq 5000 -# A cutoff in vacuum! ## electrostatics ##pme on ##pmeGridSpacing 1.2 -print "DONE [myReplica]" colvars on -colvarsConfig abf_12-15_less.colvars +colvarsConfig sabf_12-15_less.colvars +source mwABF_selection.tcl -# Selection rules. -replicaUniformPatchGrids on +# Parameters for mwABF with selection +set cycleNum 4 +set selectFreq 10000 -source selectionRules.tcl -## Selection rule parameters -set cycleNum 20 -set shareFreq 2000 -set selectCycle 5 -set repNum [numReplicas] -set rep [myReplica] -set sampleRad 0 -set binNum [cv bias abf binnum] -set percentStop 40 +run_mwABF_selection $cycleNum $selectFreq -set s "REPLICA_SETUP" -foreach var {cycleNum shareFreq selectCycle repNum rep sampleRad binNum percentStop} { - set s "$s $var [set $var]" -} -print "$s" - -for {set b 0} {$b < $binNum} {incr b} { - set binLastCount($b) [cv bias abf bincount $b] -} - -for {set i 0} {$i < $cycleNum} {incr i} { - print "CYCLE $i" - # Run the steps. - run $shareFreq - # Share the biases. - cv bias abf share - print "REPLICA_SHARED_ABF" - - if {$i % $selectCycle == 0} { - ## Get the bin. - set bin [cv bias abf bin] - # Keep within the domain. - if {$bin < 0} {set bin 0} - if {$bin >= $binNum} {set bin [expr {$binNum-1}]} - - ## Get the count in the vicinity of the bin. - set count0 [localCount $bin $binNum $sampleRad abf] - #set count [expr {$count0 - $binLastCount($bin)}] - set count $count0 -# print "REP_BIN REP $rep BIN $bin LASTCOUNT $binLastCount($bin) COUNT0 $count COUNT $count" - - replicaBarrier - if {$rep > 0} { - ## Send the count to Replica 0. - replicaSend $count 0 - ## Receive the source and destination replicas. - set srcDestList [replicaRecv 0] - } else { - ## Build the weight list. - # The weight is the inverse count. - set weightList [list [expr {1.0/($count+1)}]] - print "REPLICA_COUNT $rep $count" - set countMin $count - set countMax $count - for {set r 1} {$r < $repNum} {incr r} { - set repCount [replicaRecv $r] - print "REPLICA_COUNT $r $repCount" - if {$repCount < $countMin} { - set countMin $repCount - } - if {$repCount > $countMax} { - set countMax $repCount - } - lappend weightList [expr {1.0/$repCount}] - } - - ## Normalize the weight list. - set weightList [normalizeWeights $weightList] - set s "REPLICA_WEIGHT_LIST" - foreach w $weightList { - set s "$s [format " %.3g" $w]" - } - print $s - - print "REPLICA_MINMAX $countMin $countMax" - if {$countMin < 1} {set countMin 1} - set percentDif [expr {(100*($countMax - $countMin))/$countMax}] - - ## Generate the list of exchanges "srcDestList" - if {$percentDif < $percentStop} { - print "REPLICA_SELECTION_DISABLED $percentDif" - ## If the relative difference between the min and max counts - ##is less than the threshold, we don't do exchanges. - set srcDestList {} - } else { - print "REPLICA_SELECTION_ENABLED $percentDif" - set cloneList [resampleWalkers $weightList] - print "REPLICA_CLONE_LIST $cloneList" - - set srcDestList [resampleExchanges $cloneList] - } - - ## Replica 0 sends the srcDestList to all other replicas, - ## so they know who to receive from and who to send to. - for {set r 1} {$r < $repNum} {incr r} { - replicaSend $srcDestList $r - } - }; # End Replica 0 work - - ## Everyone should have an identical copy of the srcDestList. - if {[llength $srcDestList] > 0} { - print "REPLICA_SRC_DEST_LIST $srcDestList" - print "REPLICA_BARRIER" - replicaBarrier - - # Do the coordinate exchanges if this replica is the source or destination. - foreach srcDest $srcDestList { - set src [lindex $srcDest 0] - set dest [lindex $srcDest 1] - - if {$src == $rep} { - print "REPLICA_ATOM_SEND $dest" - replicaAtomSend $dest - } elseif {$dest == $rep} { - print "REPLICA_ATOM_RECV $src" - replicaAtomRecv $src - } - } - - # Since the resampling was successful, we set binLastCount to the current values. - for {set b 0} {$b < $binNum} {incr b} { - set binLastCount($b) [cv bias abf bincount $b] - } - } else { - print "REPLICA_SRC_DEST_EMPTY" - # The resampling did not happen. binLastCount is not updated. - } - } -} diff --git a/namd/tests/library/multiple_walker_abf/run_mw_tests.sh b/namd/tests/library/multiple_walker_abf/run_mw_tests.sh index 46ec5049f..3b365c285 100755 --- a/namd/tests/library/multiple_walker_abf/run_mw_tests.sh +++ b/namd/tests/library/multiple_walker_abf/run_mw_tests.sh @@ -11,16 +11,25 @@ if (( $# != 1 )); then fi NAMD_EXEC=$1 -#NAMD_EXEC=$HOME/Software/NAMD_2.10b1_Source_colvars/Linux-x86_64-g++/namd2 +# Runner for the netlrts version +RUNNER="${HOME}/progs/namd/Linux-x86_64-g++.netlrts/charmrun ++local +p4" + +# Runner for the MPI version +#RUNNER="mpirun -n 4" + unset l i=0 for f in mw_independent.namd mw_shared.namd mw_selection.namd; do nam=${f%.*} rm -f out_${nam}.*.log - mpirun -n 4 $NAMD_EXEC +replicas 4 $f +stdout out_${nam}.%d.log + $RUNNER $NAMD_EXEC +replicas 4 $f +stdout out_${nam}.%d.log - grad=out_${nam}.0.grad + if [ $f = "mw_independent.namd" ]; then + grad=out_${nam}.0.grad + else + grad=out_${nam}.0.all.grad + fi rmsd=$( paste -d ' ' da10_12-32_indepSaga0.0.grad $grad | awk '!/^#/ && NF==4 {n++; s+=($4-$2)*($4-$2)}; END {print sqrt(s/n)}' ) echo "RMSD $rmsd" l[$i]=$( printf "%s\t\t%.5f" $nam $rmsd ) diff --git a/namd/tests/library/multiple_walker_abf/sabf_12-15_less.colvars b/namd/tests/library/multiple_walker_abf/sabf_12-15_less.colvars index dfe3a1c3e..9b74c19e7 100644 --- a/namd/tests/library/multiple_walker_abf/sabf_12-15_less.colvars +++ b/namd/tests/library/multiple_walker_abf/sabf_12-15_less.colvars @@ -29,7 +29,6 @@ atomNumbers { 10 } } abf { - name abf colvars endToEnd fullSamples 500 historyFreq 20000 diff --git a/namd/tests/library/multiple_walker_abf/selectionRules.tcl b/namd/tests/library/multiple_walker_abf/selectionRules.tcl deleted file mode 100644 index fde136b36..000000000 --- a/namd/tests/library/multiple_walker_abf/selectionRules.tcl +++ /dev/null @@ -1,98 +0,0 @@ -# Procs for implementing selection rules. -# Author: Jeff Comer - -proc normalizeWeights {weightList} { - ## Normalize the weights. - set weightSum 0.0 - foreach weight $weightList { - set weightSum [expr {$weightSum + $weight}] - } - - set r 0 - set wnList {} - foreach weight $weightList { - lappend wnList [expr {double($weight)/$weightSum}] - incr r - } - return $wnList -} - - -# Determine the number of clones for each walker. -# Weights must be normalized. -proc resampleWalkers {weightList} { - set num [llength $weightList] - - ## Get the number of clones for each walker. - set wbar(0) [lindex $weightList 0] - set u [expr {rand()}] - set cloneList [list [expr { int($num*$wbar(0)+$u) }] ] - for {set r 1} {$r < $num} {incr r} { - set r0 [expr {$r-1}] - set wbar($r) [expr {$wbar($r0) + [lindex $weightList $r]}] - lappend cloneList [expr {int($num*$wbar($r)+$u) - int($num*$wbar($r0)+$u) }] - } - return $cloneList -} - -# Determine the minimal exchanges that must be made to resample. -proc resampleExchanges {cloneList} { - ## Make a list of exchanges. - set cloneZeroList {} - set cloneMultList {} - set r 0 - foreach cloneNum $cloneList { - if {$cloneNum == 0} { lappend cloneZeroList $r } - if {$cloneNum > 1} { lappend cloneMultList $r } - incr r - } - - # Is nothing cloned? - if {[llength $cloneZeroList] == 0} { - return {} - } - - ## Walkers cloned multiple times are copied to walkers cloned zero times. - ## Make the list of exchanges srcDestList - set srcDestList {} - set zeroInd 0 - foreach mult $cloneMultList { - # We get one clone just by leaving the walker where it is. - set extraNum [expr {[lindex $cloneList $mult]-1}] - for {set j 0} {$j < $extraNum} {incr j} { - set dest [lindex $cloneZeroList $zeroInd] - lappend srcDestList [list $mult $dest] - incr zeroInd - } - } - - return $srcDestList -} - -proc localCount {bin binNum sampleRad biasName} { - set b0 [expr {$bin-$sampleRad}] - set b1 [expr {$bin+$sampleRad}] - - set count 0 - for {set b $b0} {$b <= $b1} {incr b} { - if {$b < 0} { - set i 0 - } elseif {$b >= $binNum} { - set i [expr {$binNum-1}] - } else { - set i $b - } - incr count [cv bias $biasName bincount $i] - } - return $count -} - -# Weights are by convention normalized. -proc weightEntropy {weightList} { - set entropySum 0.0 - set num [llength $weightList] - foreach w $weightList { - set entropySum [expr {$entropySum + $w($r)*log($w($r)*$num)}] - } - return $entropySum -} diff --git a/src/colvar.cpp b/src/colvar.cpp index 0f431955d..20fdad1ac 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -1221,7 +1221,7 @@ int colvar::init_dependencies() { // Initialize feature_states for each instance feature_states.reserve(f_cv_ntot); - for (i = 0; i < f_cv_ntot; i++) { + for (i = feature_states.size(); i < f_cv_ntot; i++) { feature_states.push_back(feature_state(true, false)); // Most features are available, so we set them so // and list exceptions below @@ -1513,6 +1513,7 @@ int colvar::collect_cvc_values() cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n"); if (after_restart) { + x_old = x_restart; if (cvm::proxy->simulation_running()) { cvm::real const jump2 = dist2(x, x_restart) / (width*width); if (jump2 > 0.25) { @@ -1707,12 +1708,13 @@ int colvar::calc_colvar_properties() // Do the same if no simulation is running (eg. VMD postprocessing) if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset || !cvm::proxy->simulation_running()) { x_ext = x; + cvm::log("Initializing extended coordinate to colvar value.\n"); if (is_enabled(f_cv_reflecting_lower_boundary) && x_ext < lower_boundary) { - cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below."); + cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.\n"); x_ext = lower_boundary; } if (is_enabled(f_cv_reflecting_upper_boundary) && x_ext > upper_boundary) { - cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above."); + cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.\n"); x_ext = upper_boundary; } @@ -1722,8 +1724,18 @@ int colvar::calc_colvar_properties() // Special case of a repeated timestep (eg. multiple NAMD "run" statements) // revert values of the extended coordinate and velocity prior to latest integration if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) { - x_ext = prev_x_ext; - v_ext = prev_v_ext; + // Detect jumps due to discrete changes in coordinates (eg. in replica exchange schemes) + cvm::real const jump2 = dist2(x, x_old) / (width*width); + if (jump2 > 0.25) { + cvm::log("Detected discrete jump in colvar value from " + + cvm::to_str(x_old) + " to " + cvm::to_str(x) + ".\n"); + cvm::log("Reinitializing extended coordinate to colvar value.\n"); + x_ext = x; + } else { + cvm::log("Reinitializing extended coordinate to last value.\n"); + x_ext = prev_x_ext; + v_ext = prev_v_ext; + } } // report the restraint center as "value" // These position and velocities come from integration at the _previous timestep_ in update_forces_energy() @@ -1939,9 +1951,8 @@ int colvar::end_of_step() if (cvm::debug()) cvm::log("End of step for colvar \""+this->name+"\".\n"); - if (is_enabled(f_cv_fdiff_velocity)) { - x_old = x; - } + // Used for fdiff_velocity and for detecting jumps for extended Lagrangian colvars + x_old = x; if (is_enabled(f_cv_subtract_applied_force)) { f_old = f; diff --git a/src/colvaratoms.cpp b/src/colvaratoms.cpp index 534ca967b..977804465 100644 --- a/src/colvaratoms.cpp +++ b/src/colvaratoms.cpp @@ -286,7 +286,7 @@ int cvm::atom_group::init_dependencies() { // Initialize feature_states for each instance // default as unavailable, not enabled feature_states.reserve(f_ag_ntot); - for (i = 0; i < colvardeps::f_ag_ntot; i++) { + for (i = feature_states.size(); i < colvardeps::f_ag_ntot; i++) { feature_states.push_back(feature_state(false, false)); } diff --git a/src/colvarbias.cpp b/src/colvarbias.cpp index 2609a6469..e9a5debbf 100644 --- a/src/colvarbias.cpp +++ b/src/colvarbias.cpp @@ -167,6 +167,10 @@ int colvarbias::init_dependencies() { init_feature(f_cvb_get_total_force, "obtain_total_force", f_type_dynamic); require_feature_children(f_cvb_get_total_force, f_cv_total_force); + // Depending on back-end, we may not obtain total force at step 0 + if (!cvm::main()->proxy->total_forces_same_step()) { + exclude_feature_self(f_cvb_get_total_force, f_cvb_step_zero_data); + } init_feature(f_cvb_output_acc_work, "output_accumulated_work", f_type_user); require_feature_self(f_cvb_output_acc_work, f_cvb_apply_force); @@ -203,7 +207,7 @@ int colvarbias::init_dependencies() { // Initialize feature_states for each instance feature_states.reserve(f_cvb_ntot); - for (i = 0; i < f_cvb_ntot; i++) { + for (i = feature_states.size(); i < f_cvb_ntot; i++) { feature_states.push_back(feature_state(true, false)); // Most features are available, so we set them so // and list exceptions below @@ -437,17 +441,25 @@ int colvarbias::bin_num() cvm::error("Error: bin_num() not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; } + int colvarbias::current_bin() { cvm::error("Error: current_bin() not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; } + int colvarbias::bin_count(int /* bin_index */) { cvm::error("Error: bin_count() not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; } +int colvarbias::local_sample_count(int /* radius */) +{ + cvm::error("Error: local_sample_count() not implemented.\n"); + return COLVARS_NOT_IMPLEMENTED; +} + int colvarbias::replica_share() { cvm::error("Error: replica_share() not implemented.\n"); diff --git a/src/colvarbias.h b/src/colvarbias.h index d59746dab..1e3c6ad33 100644 --- a/src/colvarbias.h +++ b/src/colvarbias.h @@ -91,8 +91,9 @@ class colvarbias // FIXME this is currently 1D only virtual int current_bin(); //// Give the count at a given bin index. - // FIXME this is currently 1D only virtual int bin_count(int bin_index); + /// Return the average number of samples in a given "radius" around current bin + virtual int local_sample_count(int radius); //// Share information between replicas, whatever it may be. virtual int replica_share(); diff --git a/src/colvarbias_abf.cpp b/src/colvarbias_abf.cpp index 94d41e767..5c22bcdfd 100644 --- a/src/colvarbias_abf.cpp +++ b/src/colvarbias_abf.cpp @@ -14,28 +14,53 @@ #include "colvarbias_abf.h" #include "colvars_memstream.h" - colvarbias_abf::colvarbias_abf(char const *key) : colvarbias(key), b_UI_estimator(false), b_CZAR_estimator(false), pabf_freq(0), - system_force(NULL), - gradients(NULL), - samples(NULL), - pmf(NULL), - z_gradients(NULL), - z_samples(NULL), - czar_gradients(NULL), - czar_pmf(NULL), - last_gradients(NULL), - last_samples(NULL) + system_force(NULL) { - colvarproxy *proxy = cvm::main()->proxy; - if (!proxy->total_forces_same_step()) { - // Samples at step zero can not be collected - feature_states[f_cvb_step_zero_data].available = false; + // colvarproxy *proxy = cvm::main()->proxy; + init_dependencies(); +} + +int colvarbias_abf::init_dependencies() { + + // First build dependency tree for parent class (colvarbias) + colvarbias::init_dependencies(); + + int i; + if (features().size() < f_cvb_ntot) { + cvm::error("Error: dependency tree for parent class colvarbias has incorrect size " + + cvm::to_str(features().size()) + ".\n"); + return COLVARS_ERROR; + } else if (features().size() == f_cvb_ntot) { + for (i = f_cvb_ntot; i < f_cvb_abf_ntot; i++) { + modify_features().push_back(new feature); + } + + // To update the FE gradient we need to get the total force + require_feature_self(f_cvb_history_dependent, f_cvb_get_total_force); + + init_feature(f_cvb_abf_extended, "ABF on extended-Lagrangian variables", f_type_static); + + // check that everything is initialized + for (i = f_cvb_ntot; i < f_cvb_abf_ntot; i++) { + if (is_not_set(i)) { + cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description); + } + } + } + + // Initialize feature_states for each instance + for (i = f_cvb_ntot; i < f_cvb_abf_ntot; i++) { + feature_states.push_back(feature_state(true, false)); + // Most features are available, so we set them so + // and list exceptions below } + + return COLVARS_OK; } @@ -55,18 +80,11 @@ int colvarbias_abf::init(std::string const &conf) // ************* parsing general ABF options *********************** - get_keyval_feature((colvarparse *)this, conf, "applyBias", f_cvb_apply_force, true); + get_keyval_feature(this, conf, "applyBias", f_cvb_apply_force, true); if (!is_enabled(f_cvb_apply_force)){ cvm::log("WARNING: ABF biases will *not* be applied!\n"); } - get_keyval(conf, "updateBias", update_bias, true); - if (update_bias) { - enable(f_cvb_history_dependent); - } else { - cvm::log("WARNING: ABF biases will *not* be updated!\n"); - } - get_keyval(conf, "hideJacobian", hide_Jacobian, false); if (hide_Jacobian) { cvm::log("Jacobian (geometric) forces will be handled internally.\n"); @@ -74,13 +92,21 @@ int colvarbias_abf::init(std::string const &conf) cvm::log("Jacobian (geometric) forces will be included in reported free energy gradients.\n"); } - get_keyval(conf, "fullSamples", full_samples, 200); - if ( full_samples <= 1 ) full_samples = 1; - min_samples = full_samples / 2; - // full_samples - min_samples >= 1 is guaranteed + full_samples = 200; + get_keyval(conf, "fullSamples", full_samples, full_samples); + get_keyval(conf, "minSamples", min_samples, full_samples / 2); + + if (full_samples <= 1 ) { + full_samples = 1; + min_samples = 0; + } + if (min_samples >= full_samples) { + return cvm::error("Error: minSamples must be lower than fullSamples\n"); + } get_keyval(conf, "inputPrefix", input_prefix, std::vector()); + history_last_step = -1; get_keyval(conf, "historyFreq", history_freq, 0); if (history_freq != 0) { if (output_freq == 0) { @@ -93,7 +119,6 @@ int colvarbias_abf::init(std::string const &conf) } } } - b_history_files = (history_freq > 0); // shared ABF get_keyval(conf, "shared", shared_on, false); @@ -109,21 +134,22 @@ int colvarbias_abf::init(std::string const &conf) // If shared_freq is not set, we default to output_freq get_keyval(conf, "sharedFreq", shared_freq, output_freq); + if ( shared_freq && output_freq % shared_freq ) { + return cvm::error("Error: outputFreq must be a multiple of sharedFreq.\n"); + } + + // Allocate these at init time if possible + local_samples.reset(new colvar_grid_count(colvars)); + local_gradients.reset(new colvar_grid_gradient(colvars, local_samples.get())); + local_pmf.reset(new integrate_potential(colvars, local_gradients.get())); } // ************* checking the associated colvars ******************* if (num_variables() == 0) { - cvm::error("Error: no collective variables specified for the ABF bias.\n"); - return COLVARS_ERROR; + return cvm::error("Error: no collective variables specified for the ABF bias.\n"); } - if (update_bias) { - // Request calculation of total force - if(enable(f_cvb_get_total_force)) return cvm::get_error(); - } - - bool b_extended = false; size_t i; for (i = 0; i < num_variables(); i++) { @@ -135,11 +161,10 @@ int colvarbias_abf::init(std::string const &conf) colvars[i]->enable(f_cv_hide_Jacobian); } - // If any colvar is extended-system (restrained style, not external with constraint), we are running eABF - if (colvars[i]->is_enabled(f_cv_extended_Lagrangian) - && !colvars[i]->is_enabled(f_cv_external)) { - b_extended = true; - } + // If any colvar is extended-system, we need to collect the extended + // system gradient + if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)) + enable(f_cvb_abf_extended); // Cannot mix and match coarse time steps with ABF because it gives // wrong total force averages - total force needs to be averaged over @@ -155,7 +180,14 @@ int colvarbias_abf::init(std::string const &conf) // and make it just a warning if some parameter is set? } - if (b_extended) { + get_keyval(conf, "updateBias", update_bias, true); + if (update_bias) { + enable(f_cvb_history_dependent); + } else { + cvm::log("WARNING: ABF biases will *not* be updated!\n"); + } + + if (is_enabled(f_cvb_abf_extended)) { cvm::main()->cite_feature("eABF implementation"); } else { cvm::main()->cite_feature("Internal-forces free energy estimator"); @@ -185,13 +217,14 @@ int colvarbias_abf::init(std::string const &conf) cvm::log("Allocating count and free energy gradient grids.\n"); } - samples = new colvar_grid_count(colvars); - gradients = new colvar_grid_gradient(colvars); - gradients->samples = samples; - samples->has_parent_data = true; + samples.reset(new colvar_grid_count(colvars)); + gradients.reset(new colvar_grid_gradient(colvars, samples.get())); + + gradients->full_samples = full_samples; + gradients->min_samples = min_samples; - // Data for eAB F z-based estimator - if ( b_extended ) { + // Data for eABF z-based estimator + if (is_enabled(f_cvb_abf_extended)) { get_keyval(conf, "CZARestimator", b_CZAR_estimator, true); if ( b_CZAR_estimator ) { cvm::main()->cite_feature("CZAR eABF estimator"); @@ -201,13 +234,11 @@ int colvarbias_abf::init(std::string const &conf) colvarparse::parse_silent); z_bin.assign(num_variables(), 0); - z_samples = new colvar_grid_count(colvars); + z_samples.reset(new colvar_grid_count(colvars)); z_samples->request_actual_value(); - z_gradients = new colvar_grid_gradient(colvars); + z_gradients.reset(new colvar_grid_gradient(colvars, z_samples.get())); z_gradients->request_actual_value(); - z_gradients->samples = z_samples; - z_samples->has_parent_data = true; - czar_gradients = new colvar_grid_gradient(colvars); + czar_gradients.reset(new colvar_grid_gradient(colvars)); } get_keyval(conf, "integrate", b_integrate, num_variables() <= 3); // Integrate for output if d<=3 @@ -217,9 +248,9 @@ int colvarbias_abf::init(std::string const &conf) cvm::error("Error: cannot integrate free energy in dimension > 3.\n"); return COLVARS_ERROR; } - pmf = new integrate_potential(colvars, gradients); - if ( b_CZAR_estimator ) { - czar_pmf = new integrate_potential(colvars, czar_gradients); + pmf.reset(new integrate_potential(colvars, gradients.get())); + if (b_CZAR_estimator) { + czar_pmf.reset(new integrate_potential(colvars, czar_gradients.get())); } // Parameters for integrating initial (and final) gradient data get_keyval(conf, "integrateMaxIterations", integrate_iterations, 10000, colvarparse::parse_silent); @@ -230,15 +261,33 @@ int colvarbias_abf::init(std::string const &conf) get_keyval(conf, "pABFintegrateTol", pabf_integrate_tol, 1e-4, colvarparse::parse_silent); } + if (b_CZAR_estimator && shared_on && cvm::main()->proxy->replica_index() == 0) { + // The pointers below are used for outputting CZAR data + // Allocate grids for collected global data, on replica 0 only + global_z_samples = new colvar_grid_count(colvars); + global_z_gradients = new colvar_grid_gradient(colvars, global_z_samples); + global_czar_gradients = new colvar_grid_gradient(colvars); + global_czar_pmf = new integrate_potential(colvars, global_czar_gradients); + } else { + // otherwise they are just aliases for the local CZAR grids + global_z_samples = z_samples.get(); + global_z_gradients = z_gradients.get(); + global_czar_gradients = czar_gradients.get(); + global_czar_pmf = czar_pmf.get(); + } + // For shared ABF, we store a second set of grids. // This used to be only if "shared" was defined, // but now we allow calling share externally (e.g. from Tcl). - last_samples = new colvar_grid_count(colvars); - last_gradients = new colvar_grid_gradient(colvars); - last_gradients->samples = last_samples; - last_samples->has_parent_data = true; + if (b_CZAR_estimator) { + last_z_samples.reset(new colvar_grid_count(colvars)); + last_z_gradients.reset(new colvar_grid_gradient(colvars, last_z_samples.get())); + } + last_samples.reset(new colvar_grid_count(colvars)); + last_gradients.reset(new colvar_grid_gradient(colvars, last_samples.get())); shared_last_step = -1; + // If custom grids are provided, read them if ( input_prefix.size() > 0 ) { read_gradients_samples(); @@ -247,7 +296,7 @@ int colvarbias_abf::init(std::string const &conf) } // if extendedLangrangian is on, then call UI estimator - if (b_extended) { + if (is_enabled(f_cvb_abf_extended)) { get_keyval(conf, "UIestimator", b_UI_estimator, false); if (b_UI_estimator) { @@ -285,63 +334,22 @@ int colvarbias_abf::init(std::string const &conf) /// Destructor colvarbias_abf::~colvarbias_abf() { - if (samples) { - delete samples; - samples = NULL; - } - - if (gradients) { - delete gradients; - gradients = NULL; - } - - if (pmf) { - delete pmf; - pmf = NULL; - } - - if (z_samples) { - delete z_samples; - z_samples = NULL; - } - - if (z_gradients) { - delete z_gradients; - z_gradients = NULL; - } - - if (czar_gradients) { - delete czar_gradients; - czar_gradients = NULL; - } - - if (czar_pmf) { - delete czar_pmf; - czar_pmf = NULL; - } - - // shared ABF - // We used to only do this if "shared" was defined, - // but now we can call shared externally - if (last_samples) { - delete last_samples; - last_samples = NULL; - } - - if (last_gradients) { - delete last_gradients; - last_gradients = NULL; - } - - if (system_force) { - delete [] system_force; - system_force = NULL; + #define delete_if_non_null(ptr) if (ptr) { delete ptr; ptr = nullptr; } + // Only deallocate objects that are not owned by unique pointers + if (global_z_samples != z_samples.get()) { + delete_if_non_null(global_z_gradients) + delete_if_non_null(global_z_samples) + delete_if_non_null(global_czar_gradients) + delete_if_non_null(global_czar_pmf) } + if (system_force) delete[] system_force; + // Need to delete here because of additional features of the derived class + // wrt colvarbias base class + delete_features(); } /// Update the FE gradient, compute and apply biasing force -/// also output data to disk if needed int colvarbias_abf::update() { @@ -351,68 +359,81 @@ int colvarbias_abf::update() for (i = 0; i < num_variables(); i++) { bin[i] = samples->current_bin_scalar(i); } + + + // *********************************************************** + // ****** ABF Part I: update the FE gradient estimate ****** + // *********************************************************** + + if (cvm::proxy->total_forces_same_step()) { // e.g. in LAMMPS, total forces are current force_bin = bin; } - if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) { + // Share data first, so that 2d/3d PMF is refreshed using new data for mw-pABF. + // shared_on can be true with shared_freq 0 if we are sharing via script + if (shared_on && shared_freq && + shared_last_step >= 0 && // we have already collected some data + cvm::step_absolute() > shared_last_step && // time has passed since the last sharing timestep + // (avoid re-sharing at last and first ts of successive run statements) + cvm::step_absolute() % shared_freq == 0) { + // Share gradients and samples for shared ABF. + replica_share(); + } + + if (can_accumulate_data() && is_enabled(f_cvb_history_dependent)) { - if (update_bias) { -// if (b_adiabatic_reweighting) { -// // Update gradients non-locally based on conditional distribution of -// // fictitious variable TODO -// -// } else + if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) { if (samples->index_ok(force_bin)) { // Only if requested and within bounds of the grid... - for (i = 0; i < num_variables(); i++) { - // get total forces (lagging by 1 timestep) from colvars - // and subtract previous ABF force if necessary - update_system_force(i); + // get total forces (lagging by 1 timestep) from colvars + // and subtract previous ABF force if necessary + update_system_force(); + + gradients->acc_force(force_bin, system_force); + if ( b_integrate ) { + pmf->update_div_neighbors(force_bin); } - gradients->acc_force(force_bin, system_force); - if ( b_integrate ) { - pmf->update_div_neighbors(force_bin); - } } - } - if ( z_gradients && update_bias ) { - for (i = 0; i < num_variables(); i++) { - z_bin[i] = z_samples->current_bin_scalar(i); - } - if ( z_samples->index_ok(z_bin) ) { + if ( z_gradients ) { for (i = 0; i < num_variables(); i++) { - // If we are outside the range of xi, the force has not been obtained above + z_bin[i] = z_samples->current_bin_scalar(i); + } + if ( z_samples->index_ok(z_bin) ) { + // If we are outside the range of z, the force has not been obtained above // the function is just an accessor, so cheap to call again anyway - update_system_force(i); + update_system_force(); + z_gradients->acc_force(z_bin, system_force); } - z_gradients->acc_force(z_bin, system_force); } - } - if ( b_integrate ) { if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) { cvm::real err; - int iter = pmf->integrate(pabf_integrate_iterations, pabf_integrate_tol, err); - if ( iter == pabf_integrate_iterations ) { - cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(pabf_integrate_tol) - + " in " + cvm::to_str(pabf_integrate_iterations) + int iter = pmf->integrate(integrate_iterations, integrate_tol, err); + if ( iter == integrate_iterations ) { + cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol) + + " in " + cvm::to_str(integrate_iterations) + " steps. Residual error: " + cvm::to_str(err)); } - pmf->set_zero_minimum(); // TODO: do this only when necessary } } } - if (!cvm::proxy->total_forces_same_step()) { + if (!(cvm::proxy->total_forces_same_step())) { // e.g. in NAMD, total forces will be available for next timestep // hence we store the current colvar bin force_bin = bin; } + + // ****************************************************************** + // ****** ABF Part II: calculate and apply the biasing force ****** + // ****************************************************************** + + // Reset biasing forces from previous timestep for (i = 0; i < num_variables(); i++) { colvar_forces[i].reset(); @@ -421,50 +442,19 @@ int colvarbias_abf::update() // Compute and apply the new bias, if applicable if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) { - cvm::real count = cvm::real(samples->value(bin)); - cvm::real fact = 1.0; + std::vector force(num_variables()); + calc_biasing_force(force); - // Factor that ensures smooth introduction of the force - if ( count < full_samples ) { - fact = (count < min_samples) ? 0.0 : - (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); + for (size_t i = 0; i < num_variables(); i++) { + colvar_forces[i].real_value = force[i]; } + } - std::vector grad(num_variables()); - if ( pabf_freq ) { - // In projected ABF, the force is the PMF gradient estimate - pmf->vector_gradient_finite_diff(bin, grad); - } else { - // Normal ABF - gradients->vector_value(bin, grad); - } - -// if ( b_adiabatic_reweighting) { -// // Average of force according to conditional distribution of fictitious variable -// // need freshly integrated PMF, gradient TODO -// } else - if ( fact != 0.0 ) { - if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) { - // Enforce a zero-mean bias on periodic, 1D coordinates - // in other words: boundary condition is that the biasing potential is periodic - // This is enforced naturally if using integrated PMF - colvar_forces[0].real_value = fact * (grad[0] - gradients->average ()); - } else { - for (i = 0; i < num_variables(); i++) { - // subtracting the mean force (opposite of the FE gradient) means adding the gradient - colvar_forces[i].real_value = fact * grad[i]; - } - } - if (cap_force) { - for (i = 0; i < num_variables(); i++) { - if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) { - colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]); - } - } - } - } - } + // ********************************* + // ****** End of ABF proper ****** + // ********************************* + // update the output prefix; TODO: move later to setup_output() function if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) { @@ -474,10 +464,6 @@ int colvarbias_abf::update() output_prefix = cvm::output_prefix() + "." + this->name; } - if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) { - // Share gradients and samples for shared ABF. - replica_share(); - } // Prepare for the first sharing. if (shared_last_step < 0) { @@ -485,7 +471,6 @@ int colvarbias_abf::update() last_gradients->copy_grid(*gradients); last_samples->copy_grid(*samples); shared_last_step = cvm::step_absolute(); - cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".\n"); } // update UI estimator every step @@ -509,6 +494,86 @@ int colvarbias_abf::update() } + // **************************************** + // ****** Helper functions for ABF ****** + // **************************************** + + +int colvarbias_abf::update_system_force() +{ + size_t i; + // System force from atomic forces (or extended Lagrangian if applicable) + + for (i = 0; i < num_variables(); i++) { + if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) { + // this colvar is already subtracting the ABF force + system_force[i] = colvars[i]->total_force().real_value; + } else { + system_force[i] = colvars[i]->total_force().real_value + - colvar_forces[i].real_value; + } + } + return COLVARS_OK; +} + + +cvm::real colvarbias_abf::smoothing_factor(cvm::real weight) +{ + cvm::real fact = 1.0; + if ( weight < full_samples ) { + if ( weight < min_samples) { + fact = 0.0; + } else { + fact = (weight - min_samples) / cvm::real(full_samples - min_samples); + } + } + return fact; +} + + +int colvarbias_abf::calc_biasing_force(std::vector &force) +{ + size_t i; + + // Pick between different types of biasing force + if ( pabf_freq ) { + // In projected ABF, the force is the PMF gradient estimate + pmf->vector_gradient_finite_diff(bin, force); + // Calculate ramp factor that ensures smooth introduction of the force + const cvm::real count = samples->value(bin); + const cvm::real fact = smoothing_factor(count); + for (i = 0; i < num_variables(); i++) { + force[i] *= fact; + } + } else { + // Normal ABF or eABF: use accumulated gradient average + gradients->vector_value_smoothed(bin, &force[0], true); + if ( (num_variables() == 1) && gradients->periodic[0] ) { + // Enforce a zero-mean bias on periodic, 1D coordinates + // in other words: boundary condition is that the biasing potential is periodic + // Only plain ABF needs this + force[0] = force[0] - gradients->average(); + } + } + + if (cap_force) { + for (i = 0; i < num_variables(); i++) { + if ( force[i] * force[i] > max_force[i] * max_force[i] ) { + force[i] = (force[i] > 0 ? max_force[i] : -1.0 * max_force[i]); + } + } + } + + return COLVARS_OK; +} + + + // ************************************ + // ****** Shared ABF functions ****** + // ************************************ + + + int colvarbias_abf::replica_share() { colvarproxy *proxy = cvm::main()->proxy; @@ -522,53 +587,82 @@ int colvarbias_abf::replica_share() { cvm::error("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n"); return COLVARS_ERROR; } + shared_on = true; // If called by a script, inform the rest of the code that we're sharing, eg. CZAR // Share gradients for shared ABF. cvm::log("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) ); + if (!local_samples) { + // We arrive here if sharing has just been enabled by a script + // in which case local arrays have not been initialized yet + local_samples.reset(new colvar_grid_count(colvars)); + local_gradients.reset(new colvar_grid_gradient(colvars, local_samples.get())); + local_pmf.reset(new integrate_potential(colvars, local_gradients.get())); + } + // Calculate the delta gradient and count for the local replica + last_gradients->delta_grid(*gradients); + // Add the delta gradient and count to the accumulated local data + local_gradients->add_grid(*last_gradients); + + last_samples->delta_grid(*samples); + local_samples->add_grid(*last_samples); + + // Count of data items. - size_t data_n = gradients->raw_data_num(); - size_t samp_start = data_n*sizeof(cvm::real); - size_t msg_total = data_n*sizeof(size_t) + samp_start; + size_t samples_n = samples->raw_data_num(); + size_t gradients_n = gradients->raw_data_num(); + + size_t samp_start = gradients_n*sizeof(cvm::real); + int msg_total = samples_n * sizeof(size_t) + samp_start; char* msg_data = new char[msg_total]; - if (proxy->replica_index() == 0) { + if (cvm::main()->proxy->replica_index() == 0) { int p; // Replica 0 collects the delta gradient and count from the others. for (p = 1; p < proxy->num_replicas(); p++) { // Receive the deltas. - proxy->replica_comm_recv(msg_data, msg_total, p); + if (proxy->replica_comm_recv(msg_data, msg_total, p) != msg_total) { + cvm::error("Error getting shared ABF data from replica."); + return COLVARS_ERROR; + } // Map the deltas from the others into the grids. + // Re-use last_gradients as temp array, erasing its contents each time last_gradients->raw_data_in((cvm::real*)(&msg_data[0])); - last_samples->raw_data_in((size_t*)(&msg_data[samp_start])); - // Combine the delta gradient and count of the other replicas // with Replica 0's current state (including its delta). gradients->add_grid( *last_gradients ); - samples->add_grid( *last_samples ); + + last_samples->raw_data_in((size_t*)(&msg_data[samp_start])); + samples->add_grid(*last_samples); } // Now we must send the combined gradient to the other replicas. gradients->raw_data_out((cvm::real*)(&msg_data[0])); samples->raw_data_out((size_t*)(&msg_data[samp_start])); + for (p = 1; p < proxy->num_replicas(); p++) { - proxy->replica_comm_send(msg_data, msg_total, p); + if (proxy->replica_comm_send(msg_data, msg_total, p) != msg_total) { + cvm::error("Error sending shared ABF data to replica."); + return COLVARS_ERROR; + } } } else { // All other replicas send their delta gradient and count. - // Calculate the delta gradient and count. - last_gradients->delta_grid(*gradients); - last_samples->delta_grid(*samples); - // Cast the raw char data to the gradient and samples. last_gradients->raw_data_out((cvm::real*)(&msg_data[0])); last_samples->raw_data_out((size_t*)(&msg_data[samp_start])); - proxy->replica_comm_send(msg_data, msg_total, 0); + if (proxy->replica_comm_send(msg_data, msg_total, 0) != msg_total) { + cvm::error("Error sending shared ABF data to replica."); + return COLVARS_ERROR; + } // We now receive the combined gradient from Replica 0. - proxy->replica_comm_recv(msg_data, msg_total, 0); + if (proxy->replica_comm_recv(msg_data, msg_total, 0) != msg_total) { + cvm::error("Error getting shared ABF data from replica 0."); + return COLVARS_ERROR; + } // We sync to the combined gradient computed by Replica 0. gradients->raw_data_in((cvm::real*)(&msg_data[0])); samples->raw_data_in((size_t*)(&msg_data[samp_start])); @@ -585,14 +679,103 @@ int colvarbias_abf::replica_share() { last_samples->copy_grid(*samples); shared_last_step = cvm::step_absolute(); + cvm::log("RMSD btw. local and global ABF gradients: " + cvm::to_str(gradients->grid_rmsd(*local_gradients))); + if (b_integrate) { - // Update divergence to account for newly shared gradients + cvm::real err; + + // Update whole divergence field to account for newly shared gradients pmf->set_div(); + pmf->integrate(integrate_iterations, integrate_tol, err); + pmf->set_zero_minimum(); + local_pmf->set_div(); + local_pmf->integrate(integrate_iterations, integrate_tol, err); + local_pmf->set_zero_minimum(); + cvm::log("RMSD btw. local and global ABF FES: " + cvm::to_str(pmf->grid_rmsd(*local_pmf))); + } + return COLVARS_OK; +} + + +int colvarbias_abf::replica_share_CZAR() { + colvarproxy *proxy = cvm::main()->proxy; + + cvm::log("shared eABF: Gathering CZAR gradient and samples from replicas at step "+cvm::to_str(cvm::step_absolute()) ); + + // Calculate the delta gradient and count for the local replica + last_z_gradients->delta_grid(*z_gradients); + last_z_samples->delta_grid(*z_samples); + + // Count of data items. + size_t samples_n = z_samples->raw_data_num(); + size_t gradients_n = z_gradients->raw_data_num(); + + size_t samp_start = gradients_n*sizeof(cvm::real); + int msg_total = samples_n*sizeof(size_t) + samp_start; + char* msg_data = new char[msg_total]; + + if (cvm::main()->proxy->replica_index() == 0) { + if (!global_z_samples) { + // We arrive here if sharing has just been enabled by a script + // Allocate grids for collective data, on replica 0 only + // overriding CZAR grids that are equal to local ones by default + global_z_samples = new colvar_grid_count(colvars); + global_z_gradients = new colvar_grid_gradient(colvars, global_z_samples); + global_czar_gradients = new colvar_grid_gradient(colvars); + global_czar_pmf = new integrate_potential(colvars, global_czar_gradients); + } + // Start by adding new data for this cycle from replica 0 + global_z_gradients->add_grid( *last_z_gradients ); + global_z_samples->add_grid( *last_z_samples ); + int p; + // Replica 0 collects the delta gradient and count from the others. + for (p = 1; p < proxy->num_replicas(); p++) { + // Receive the deltas. + if (proxy->replica_comm_recv(msg_data, msg_total, p) != msg_total) { + cvm::error("Error getting shared ABF data from replica."); + return COLVARS_ERROR; + } + + // Map the deltas from the others into the grids. + // Re-use last_z_gradients, erasing its contents each time + last_z_gradients->raw_data_in((cvm::real*)(&msg_data[0])); + last_z_samples->raw_data_in((size_t*)(&msg_data[samp_start])); + + // Combine the new gradient and count of the other replicas + // with Replica 0's current state (including its delta). + global_z_gradients->add_grid( *last_z_gradients ); + global_z_samples->add_grid( *last_z_samples ); + } + + } else { + // All other replicas send their delta gradient and count. + last_z_gradients->raw_data_out((cvm::real*)(&msg_data[0])); + last_z_samples->raw_data_out((size_t*)(&msg_data[samp_start])); + if (proxy->replica_comm_send(msg_data, msg_total, 0) != msg_total) { + cvm::error("Error sending shared ABF data to replica."); + return COLVARS_ERROR; + } } + + // Without a barrier it's possible that one replica starts + // share 2 when other replicas haven't finished share 1. + proxy->replica_comm_barrier(); + // Done syncing the replicas. + delete[] msg_data; + + // Copy the current gradient and count values into last. + last_z_gradients->copy_grid(*z_gradients); + last_z_samples->copy_grid(*z_samples); + return COLVARS_OK; } + // ***************************** + // ****** I/O functions ****** + // ***************************** + + size_t colvarbias_abf::replica_share_freq() const { return shared_freq; @@ -638,45 +821,77 @@ template int colvarbias_abf::write_grid_to_file(T const *grid, } -void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close) +void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close, bool local) { colvarproxy *proxy = cvm::main()->proxy; - write_grid_to_file(samples, prefix + ".count", close); - write_grid_to_file(gradients, prefix + ".grad", close); + // The following are local aliases for the class' unique pointers + colvar_grid_count *samples_out, *z_samples_out; + colvar_grid_gradient *gradients_out, *z_gradients_out, *czar_gradients_out; + integrate_potential *pmf_out, *czar_pmf_out; + + // In shared ABF, write grids containing local data only if requested + if (local) { + samples_out = local_samples.get(); + gradients_out = local_gradients.get(); + pmf_out = local_pmf.get(); + // Update the divergence before integrating the local PMF below + // only needs to happen here, just before output + local_pmf->set_div(); + z_samples_out = z_samples.get(); + z_gradients_out = z_gradients.get(); + czar_gradients_out = czar_gradients.get(); + czar_pmf_out = czar_pmf.get(); + } else { + samples_out = samples.get(); + gradients_out = gradients.get(); + pmf_out = pmf.get(); + // Note: outside of shared ABF, "global" CZAR grids are just the local ones + z_samples_out = global_z_samples; + z_gradients_out = global_z_gradients; + czar_gradients_out = global_czar_gradients; + czar_pmf_out = global_czar_pmf; + } + + write_grid_to_file(samples_out, prefix + ".count", close); + write_grid_to_file(gradients_out, prefix + ".grad", close); if (b_integrate) { // Do numerical integration (to high precision) and output a PMF cvm::real err; - pmf->integrate(integrate_iterations, integrate_tol, err); - pmf->set_zero_minimum(); - write_grid_to_file(pmf, prefix + ".pmf", close); + // Divergence has already been updated on the fly for the global PMF (member data 'pmf') + pmf_out->integrate(integrate_iterations, integrate_tol, err); + pmf_out->set_zero_minimum(); + write_grid_to_file(pmf_out, prefix + ".pmf", close); } if (b_CZAR_estimator) { // Write eABF CZAR-related quantities - write_grid_to_file(z_samples, prefix + ".zcount", close); + write_grid_to_file(z_samples_out, prefix + ".zcount", close); if (b_czar_window_file) { - write_grid_to_file(z_gradients, prefix + ".zgrad", close); + write_grid_to_file(z_gradients_out, prefix + ".zgrad", close); } - // Calculate CZAR estimator of gradients - for (std::vector ix = czar_gradients->new_index(); - czar_gradients->index_ok(ix); czar_gradients->incr(ix)) { - for (size_t n = 0; n < czar_gradients->multiplicity(); n++) { - czar_gradients->set_value(ix, z_gradients->value_output(ix, n) - - proxy->target_temperature() * proxy->boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n); + // Update the CZAR estimator of gradients, except at step 0 + // in which case we preserve any existing data (e.g. read via inputPrefix, used to join strata in stratified eABF) + if (cvm::step_relative() > 0) { + for (std::vector iz_bin = czar_gradients_out->new_index(); + czar_gradients_out->index_ok(iz_bin); czar_gradients_out->incr(iz_bin)) { + for (size_t n = 0; n < czar_gradients_out->multiplicity(); n++) { + czar_gradients_out->set_value(iz_bin, z_gradients_out->value_output(iz_bin, n) + - proxy->target_temperature() * proxy->boltzmann() * z_samples_out->log_gradient_finite_diff(iz_bin, n), n); + } } } - write_grid_to_file(czar_gradients, prefix + ".czar.grad", close); + write_grid_to_file(czar_gradients_out, prefix + ".czar.grad", close); if (b_integrate) { // Do numerical integration (to high precision) and output a PMF cvm::real err; - czar_pmf->set_div(); - czar_pmf->integrate(integrate_iterations, integrate_tol, err); - czar_pmf->set_zero_minimum(); - write_grid_to_file(czar_pmf, prefix + ".czar.pmf", close); + czar_pmf_out->set_div(); + czar_pmf_out->integrate(integrate_iterations, integrate_tol, err); + czar_pmf_out->set_zero_minimum(); + write_grid_to_file(czar_pmf_out, prefix + ".czar.pmf", close); } } return; @@ -686,59 +901,67 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool clo // For Tcl implementation of selection rules. /// Give the total number of bins for a given bias. int colvarbias_abf::bin_num() { - return samples->number_of_points(0); + return samples->number_of_points(); } + /// Calculate the bin index for a given bias. int colvarbias_abf::current_bin() { - return samples->current_bin_scalar(0); + return samples->current_bin_flat_bound(); } + /// Give the count at a given bin index. int colvarbias_abf::bin_count(int bin_index) { if (bin_index < 0 || bin_index >= bin_num()) { cvm::error("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index)); return -1; } - std::vector ix(1,(int)bin_index); - return samples->value(ix); + return int(samples->get_value(bin_index)); } +// Return the average number of samples in a given "radius" around current bin +int colvarbias_abf::colvarbias_abf::local_sample_count(int radius) { + return samples->local_sample_count(radius); +} int colvarbias_abf::read_gradients_samples() { - int error_code = COLVARS_OK; + int err = COLVARS_OK; + // Reading the CZAR gradients is necessary for joining strata in stratified eABF + std::unique_ptr czar_gradients_in; - std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name; + if (b_CZAR_estimator) { + // CZAR gradients are usually computed as needed from z-gradients and z_samples + // Therefore the czar_gradients grid is not linked to a sampling grid + // Here we define a temporary czar_gradients grid linked to z_samples, + // to correctly average input gradients if overlapping + czar_gradients_in.reset(new colvar_grid_gradient(colvars, z_samples.get())); + } for ( size_t i = 0; i < input_prefix.size(); i++ ) { - samples_in_name = input_prefix[i] + ".count"; - gradients_in_name = input_prefix[i] + ".grad"; - z_samples_in_name = input_prefix[i] + ".zcount"; - z_gradients_in_name = input_prefix[i] + ".zgrad"; - // For user-provided files, the per-bias naming scheme may not apply - cvm::log("Reading sample count from " + samples_in_name + - " and gradient from " + gradients_in_name); - - error_code |= samples->read_multicol(samples_in_name, - "ABF samples file", - true); + std::string prefix = input_prefix[i]; - error_code |= gradients->read_multicol(gradients_in_name, - "ABF gradient file", - true); + // For user-provided files, the per-bias naming scheme may not apply + err |= samples->read_multicol(prefix + ".count", "ABF samples file", true); + err |= gradients->read_multicol(prefix + ".grad", "ABF gradient file", true); if (b_CZAR_estimator) { // Read eABF z-averaged data for CZAR - cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name); - error_code |= z_samples->read_multicol(z_samples_in_name, - "eABF z-histogram file", - true); - error_code |= z_gradients->read_multicol(z_gradients_in_name, - "eABF z-gradient file", - true); + err |= z_samples->read_multicol(prefix + ".zcount", "eABF z-histogram file", true); + err |= z_gradients->read_multicol(prefix + ".zgrad", "eABF z-gradient file", true); + err |= czar_gradients_in->read_multicol(prefix + ".czar.grad", "eABF CZAR gradient file", true); } } - return error_code; + if (b_CZAR_estimator) { + // Now copy real CZAR gradients (divided by total count) to the final grid + for (std::vector ix = czar_gradients->new_index(); + czar_gradients->index_ok(ix); czar_gradients->incr(ix)) { + for (size_t n = 0; n < czar_gradients->multiplicity(); n++) { + czar_gradients->set_value(ix, czar_gradients_in->value_output(ix, n), n); + } + } + } + return err; } @@ -747,9 +970,9 @@ template OST & colvarbias_abf::write_state_data_template_(OST &os auto flags = os.flags(); os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format + write_state_data_key(os, "samples"); samples->write_raw(os, 8); - os.flags(flags); write_state_data_key(os, "gradient"); gradients->write_raw(os, 8); @@ -758,7 +981,6 @@ template OST & colvarbias_abf::write_state_data_template_(OST &os os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format write_state_data_key(os, "z_samples"); z_samples->write_raw(os, 8); - os.flags(flags); write_state_data_key(os, "z_gradient"); z_gradients->write_raw(os, 8); } @@ -843,19 +1065,34 @@ int colvarbias_abf::write_output_files() cvm::log("ABF bias trying to write gradients and samples to disk"); } - if (shared_on && cvm::main()->proxy->replica_index() > 0 - && ! (b_CZAR_estimator || b_UI_estimator) ) { - // No need to report the same data as replica 0, let it do the I/O job - // except if using an eABF FE estimator - return COLVARS_OK; + // In shared eABF/CZAR, the communication routine needs to run on all ranks + if (shared_on) { + if (b_CZAR_estimator) replica_share_CZAR(); } - write_gradients_samples(output_prefix); - if (b_history_files) { - if ((cvm::step_absolute() % history_freq) == 0) { - write_gradients_samples(output_prefix + ".hist", false); + // In shared setting, output local data for all replicas + if (shared_on) { + // Write local data on all replicas + write_gradients_samples(output_prefix, true, true); + + if (cvm::main()->proxy->replica_index() > 0) { + // No need to report the same data as replica 0, let it do the I/O job + return COLVARS_OK; } } + // In shared ABF, only replica 0 reaches this + // filename prefix for master replica + // used in mwABF to distinguish local from complete data + std::string master_prefix = (shared_on ? output_prefix + ".all" : output_prefix); + write_gradients_samples(master_prefix); + + if ((history_freq > 0) && + (!shared_on || cvm::main()->proxy->replica_index() == 0) && // if shared, only on replica 0 + (cvm::step_absolute() % history_freq == 0) && // at requested frequency + (cvm::step_absolute() != history_last_step)) { // not twice the same timestep + write_gradients_samples(master_prefix + ".hist", false); + history_last_step = cvm::step_absolute(); + } if (b_UI_estimator) { eabf_UI.calc_pmf(); @@ -873,14 +1110,13 @@ int colvarbias_abf::calc_energy(std::vector const *values) if (num_variables() > 1 || values != NULL) { // Use simple estimate: neglect effect of fullSamples, // return value at center of bin - if (pmf != NULL) { - std::vector const curr_bin = values ? + if (pmf) { + std::vector curr_bin = values ? pmf->get_colvars_index(*values) : pmf->get_colvars_index(); - - if (pmf->index_ok(curr_bin)) { - bias_energy = pmf->value(curr_bin); - } + pmf->set_zero_minimum(); + pmf->wrap_to_edge(curr_bin, curr_bin); // Closest edge if outside of grid + bias_energy = pmf->value(curr_bin); } return COLVARS_OK; } @@ -895,30 +1131,21 @@ int colvarbias_abf::calc_energy(std::vector const *values) cvm::real sum = 0.0; for (int i = 0; i < home; i++) { std::vector ix(1,i); - - // Include the full_samples factor if necessary. - unsigned int count = samples->value(ix); - cvm::real fact = 1.0; - if ( count < full_samples ) { - fact = (count < min_samples) ? 0.0 : - (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); - } - if (count > 0) sum += fact*gradients->value(ix)/count*gradients->widths[0]; + // Include the smoothing factor if necessary. + sum += gradients->value_output_smoothed(ix, true) * gradients->widths[0]; } // Integrate the gradient up to the current position in the home interval, a fractional portion of a bin. std::vector ix(1,home); cvm::real frac = gradients->current_bin_scalar_fraction(0); - unsigned int count = samples->value(ix); - cvm::real fact = 1.0; - if ( count < full_samples ) { - fact = (count < min_samples) ? 0.0 : - (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); - } - if (count > 0) - sum += fact*gradients->value(ix)/count*gradients->widths[0]*frac; + sum += gradients->value_output_smoothed(ix, true) * gradients->widths[0] * frac; // The applied potential is the negative integral of force samples. bias_energy = -sum; return COLVARS_OK; } + + +// Static members + +std::vector colvarbias_abf::cvb_abf_features; diff --git a/src/colvarbias_abf.h b/src/colvarbias_abf.h index ede514d6e..ac645146f 100644 --- a/src/colvarbias_abf.h +++ b/src/colvarbias_abf.h @@ -14,13 +14,14 @@ #include #include #include +#include #include "colvarproxy.h" #include "colvarbias.h" #include "colvargrid.h" #include "colvar_UIestimator.h" -typedef cvm::real* gradient_t; +typedef cvm::real *gradient_t; /// ABF bias @@ -54,8 +55,8 @@ class colvarbias_abf : public colvarbias { size_t full_samples; /// Number of samples per bin before applying a scaled-down biasing force size_t min_samples; - /// Write combined files with a history of all output data? - bool b_history_files; + /// Latest absolute time step at which history files were written + cvm::step_number history_last_step; /// Write CZAR output file for stratified eABF (.zgrad) bool b_czar_window_file; /// Number of timesteps between recording data in history files (if non-zero) @@ -96,38 +97,28 @@ class colvarbias_abf : public colvarbias { gradient_t system_force; /// n-dim grid of free energy gradients - colvar_grid_gradient *gradients; + std::unique_ptr gradients; /// n-dim grid of number of samples - colvar_grid_count *samples; + std::unique_ptr samples; /// n-dim grid of pmf (dimension 1 to 3) - integrate_potential *pmf; + std::unique_ptr pmf; /// n-dim grid: average force on "real" coordinate for eABF z-based estimator - colvar_grid_gradient *z_gradients; + std::unique_ptr z_gradients; /// n-dim grid of number of samples on "real" coordinate for eABF z-based estimator - colvar_grid_count *z_samples; - /// n-dim grid containing CZAR estimator of "real" free energy gradients - colvar_grid_gradient *czar_gradients; + std::unique_ptr z_samples; + /// n-dim grid containing CZAR estimatr of "real" free energy gradients + std::unique_ptr czar_gradients; /// n-dim grid of CZAR pmf (dimension 1 to 3) - integrate_potential *czar_pmf; + std::unique_ptr czar_pmf; - inline int update_system_force(size_t i) - { - if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) { - // this colvar is already subtracting the ABF force - system_force[i] = colvars[i]->total_force().real_value; - } else { - system_force[i] = colvars[i]->total_force().real_value - - colvar_forces[i].real_value; - // If hideJacobian is active then total_force has an extra term of -fj - // which is the Jacobian-compensating force at the colvar level - } - if (cvm::debug()) - cvm::log("ABF System force calc: cv " + cvm::to_str(i) + - " fs " + cvm::to_str(system_force[i]) + - " = ft " + cvm::to_str(colvars[i]->total_force().real_value) + - " - fa " + cvm::to_str(colvar_forces[i].real_value)); - return COLVARS_OK; - } + /// Calculate system force for all colvars + int update_system_force(); + + /// Calulate the biasing force for the current bin + int calc_biasing_force(std::vector &force); + + /// Calulate the smoothing factor to apply to biasing forces for given local count + cvm::real smoothing_factor(cvm::real weight); // shared ABF bool shared_on; @@ -137,11 +128,29 @@ class colvarbias_abf : public colvarbias { // Share between replicas -- may be called independently of update virtual int replica_share(); + // Share data needed for CZAR between replicas - called before output only + int replica_share_CZAR(); + + /// Report the frequency at which this bias needs to communicate with replicas virtual size_t replica_share_freq() const; - // Store the last set for shared ABF - colvar_grid_gradient *last_gradients; - colvar_grid_count *last_samples; + // Data just after the last share (start of cycle) in shared ABF + std::unique_ptr last_gradients; + std::unique_ptr last_samples; + // eABF/CZAR local data last shared + std::unique_ptr last_z_gradients; + std::unique_ptr last_z_samples; + // ABF data from local replica only in shared ABF + std::unique_ptr local_gradients; + std::unique_ptr local_samples; + std::unique_ptr local_pmf; + // eABF/CZAR data collected from all replicas in shared eABF on replica 0 + // if non-shared, aliases of regular CZAR grids, for output purposes + colvar_grid_gradient *global_z_gradients; + colvar_grid_count *global_z_samples; + colvar_grid_gradient *global_czar_gradients; + integrate_potential *global_czar_pmf; + // For Tcl implementation of selection rules. /// Give the total number of bins for a given bias. @@ -150,14 +159,17 @@ class colvarbias_abf : public colvarbias { virtual int current_bin(); //// Give the count at a given bin index. virtual int bin_count(int bin_index); + /// Return the average number of samples in a given "radius" around current bin + virtual int local_sample_count(int radius); /// Write human-readable FE gradients and sample count, and DX file in dim > 2 - void write_gradients_samples(const std::string &prefix, bool close = true); + /// \param local write grids contining replica-local data in shared ABF + void write_gradients_samples(const std::string &prefix, bool close = true, bool local = false); /// Read human-readable FE gradients and sample count (if not using restart) int read_gradients_samples(); - /// Template used in write_gradient_samples() + /// Shorthand template used in write_gradient_samples() template int write_grid_to_file(T const *grid, std::string const &name, bool close); @@ -184,6 +196,38 @@ class colvarbias_abf : public colvarbias { /// Calculate the bias energy for 1D ABF virtual int calc_energy(std::vector const *values); -}; + /// Initialize specific dependencies of ABF derived class + /// Adding them to those of base class colvarbias + /// Alternately we could overload the init_dependencies() function + virtual int init_dependencies(); + + enum features_bias_abf { + /// Start after generic cvb features + /// There is at least one ext-Lagrangian colvar -> run eABF + f_cvb_abf_extended = f_cvb_ntot, + /// Total number of features for an ABF bias + f_cvb_abf_ntot + }; + + /// \brief Implementation of the feature list for colvarbias_abf + static std::vector cvb_abf_features; + + /// \brief Implementation of the feature list accessor for colvarbias + virtual const std::vector &features() const + { + return cvb_abf_features; + } + virtual std::vector &modify_features() + { + return cvb_abf_features; + } + static void delete_features() { + for (size_t i=0; i < cvb_abf_features.size(); i++) { + delete cvb_abf_features[i]; + } + cvb_abf_features.clear(); + } + +}; #endif diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 0d05e176a..4c7ab09dc 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -303,7 +303,7 @@ int colvar::cvc::init_dependencies() { // default as available, not enabled // except dynamic features which default as unavailable feature_states.reserve(f_cvc_ntot); - for (i = 0; i < colvardeps::f_cvc_ntot; i++) { + for (i = feature_states.size(); i < colvardeps::f_cvc_ntot; i++) { bool avail = is_dynamic(i) ? false : true; feature_states.push_back(feature_state(avail, false)); } diff --git a/src/colvardeps.cpp b/src/colvardeps.cpp index 2a22b1cb5..46b791756 100644 --- a/src/colvardeps.cpp +++ b/src/colvardeps.cpp @@ -129,6 +129,11 @@ int colvardeps::enable(int feature_id, int res; size_t i, j; bool ok; + + if (feature_id < 0 || feature_id >= int(features().size())) { + cvm::error("Error: colvardeps::enable() called with invalid feature_id " + cvm::to_str(feature_id) + "\n"); + return COLVARS_ERROR; + } feature *f = features()[feature_id]; feature_state *fs = &feature_states[feature_id]; diff --git a/src/colvargrid.cpp b/src/colvargrid.cpp index 5aefbb8c4..505ed286b 100644 --- a/src/colvargrid.cpp +++ b/src/colvargrid.cpp @@ -297,30 +297,63 @@ cvm::real colvar_grid_scalar::entropy() const return bin_volume * sum; } +/// \brief Return the RMSD between this grid's data and another one +/// Grids must have the same dimensions. +cvm::real colvar_grid_scalar::grid_rmsd(colvar_grid_scalar const &other_grid) const +{ + if (other_grid.data.size() != this->data.size()) { + cvm::error("Error: trying to subtract two grids with " + "different size.\n"); + return -1.; + } + + cvm::real sum2 = 0.0; + + if (samples && other_grid.samples) { + for (size_t i = 0; i < data.size(); i++) { + size_t n = samples->get_value(i); + cvm::real us = n ? data[i] / n : 0.0; + n = other_grid.samples->get_value(i); + cvm::real them = n ? other_grid.data[i ] / n : 0.0; + cvm::real d = us - them; + sum2 += d*d; + } + } else { + for (size_t i = 0; i < data.size(); i++) { + cvm::real d = other_grid.data[i] - data[i]; + sum2 += d*d; + } + } + + return sqrt(sum2/this->data.size()); +} + colvar_grid_gradient::colvar_grid_gradient() - : colvar_grid(), - samples(NULL), - weights(NULL) + : colvar_grid(), samples(NULL), full_samples(0), min_samples(0) {} + colvar_grid_gradient::colvar_grid_gradient(std::vector const &nx_i) - : colvar_grid(nx_i, 0.0, nx_i.size()), - samples(NULL), - weights(NULL) + : colvar_grid(nx_i, 0.0, nx_i.size()), samples(NULL), full_samples(0), min_samples(0) {} + colvar_grid_gradient::colvar_grid_gradient(std::vector &colvars) - : colvar_grid(colvars, 0.0, colvars.size()), - samples(NULL), - weights(NULL) + : colvar_grid(colvars, 0.0, colvars.size()), samples(NULL), full_samples(0), min_samples(0) {} +colvar_grid_gradient::colvar_grid_gradient(std::vector &colvars, colvar_grid_count *samples_in) + : colvar_grid(colvars, 0.0, colvars.size()), samples(samples_in), full_samples(0), min_samples(0) +{ + samples_in->has_parent_data = true; +} + + colvar_grid_gradient::colvar_grid_gradient(std::string &filename) : colvar_grid(), - samples(NULL), - weights(NULL) + samples(NULL) { std::istream &is = cvm::main()->proxy->input_stream(filename, "gradient file"); @@ -524,9 +557,38 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os) } +/// \brief Return the RMSD between this grid's data and another one +/// Grids must have the same dimensions. +cvm::real colvar_grid_gradient::grid_rmsd(colvar_grid_gradient const &other_grid) const +{ + if (other_grid.multiplicity() != this->multiplicity()) { + cvm::error("Error: trying to subtract two grids with " + "different multiplicity.\n"); + return -1.; + } + + if (other_grid.data.size() != this->data.size()) { + cvm::error("Error: trying to subtract two grids with " + "different size.\n"); + return -1.; + } + + cvm::real sum2 = 0.0; + std::vector ix; + size_t imult; + for (ix = new_index(); index_ok(ix); incr(ix)) { + for (imult = 0; imult < this->multiplicity(); imult++) { + cvm::real d = this->value_output(ix, imult) - other_grid.value_output(ix, imult); + sum2 += d*d; + } + } + return sqrt(sum2/this->data.size()); +} + integrate_potential::integrate_potential(std::vector &colvars, colvar_grid_gradient * gradients) : colvar_grid_scalar(colvars, true), + b_smoothed(false), gradients(gradients) { // parent class colvar_grid_scalar is constructed with margin option set to true @@ -556,7 +618,8 @@ integrate_potential::integrate_potential(std::vector &colvars, colvar_ integrate_potential::integrate_potential(colvar_grid_gradient * gradients) - : gradients(gradients) + : b_smoothed(false), + gradients(gradients) { nd = gradients->num_variables(); nx = gradients->number_of_points_vec(); @@ -578,7 +641,7 @@ integrate_potential::integrate_potential(colvar_grid_gradient * gradients) } -int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::real & err) +int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::real & err, bool verbose) { int iter = 0; @@ -591,22 +654,24 @@ int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::r } else { corr = 0.0; } - std::vector ix; // Iterate over valid indices in gradient grid for (ix = new_index(); gradients->index_ok(ix); incr(ix)) { set_value(ix, sum); - sum += (gradients->value_output(ix) - corr) * widths[0]; + cvm::real val = gradients->value_output_smoothed(ix, b_smoothed); + sum += (val - corr) * widths[0]; } if (index_ok(ix)) { // This will happen if non-periodic: then PMF grid has one extra bin wrt gradient grid + // If not, sum should be zero set_value(ix, sum); } } else if (nd <= 3) { nr_linbcg_sym(divergence, data, tol, itmax, iter, err); - cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err) + "\n"); + if (verbose) + cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err)); } else { cvm::error("Cannot integrate PMF in dimension > 3\n"); @@ -630,7 +695,7 @@ void integrate_potential::update_div_neighbors(const std::vector &ix0) std::vector ix(ix0); int i, j, k; - // If not periodic, expanded grid ensures that neighbors of ix0 are valid grid points + // If not periodic, expanded grid ensures that upper neighbors of ix0 are valid grid points if (nd == 1) { return; @@ -662,29 +727,23 @@ void integrate_potential::update_div_neighbors(const std::vector &ix0) } } + void integrate_potential::get_grad(cvm::real * g, std::vector &ix) { - size_t count, i; - bool edge = gradients->wrap_edge(ix); - if (!edge) { - if (gradients->samples) { - count = gradients->samples->value(ix); - } else { - count = 1; - } - cvm::real fact = 0.0; - if (count) fact = 1.0 / count; - cvm::real const *grad = &(gradients->value(ix)); + size_t i; + bool edge = gradients->wrap_detect_edge(ix); // Detect edge if non-PBC - for ( i = 0; ivector_value_smoothed(ix, g, b_smoothed); } + void integrate_potential::update_div_local(const std::vector &ix0) { const size_t linear_index = address(ix0); diff --git a/src/colvargrid.h b/src/colvargrid.h index ef3d96734..94a1a2534 100644 --- a/src/colvargrid.h +++ b/src/colvargrid.h @@ -374,7 +374,7 @@ template class colvar_grid : public colvarparse { { for (size_t i = 0; i < nd; i++) { if (periodic[i]) { - ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result + ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined) } else { if (ix[i] < 0 || ix[i] >= nx[i]) { cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: " @@ -387,25 +387,58 @@ template class colvar_grid : public colvarparse { /// Wrap an index vector around periodic boundary conditions /// or detects edges if non-periodic - inline bool wrap_edge(std::vector & ix) const + inline bool wrap_detect_edge(std::vector & ix) const { bool edge = false; for (size_t i = 0; i < nd; i++) { if (periodic[i]) { - ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result - } else if (ix[i] == -1 || ix[i] == nx[i]) { + ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined) + } else if (ix[i] < 0 || ix[i] >= nx[i]) { edge = true; } } return edge; } + /// Wrap an index vector around periodic boundary conditions + /// or brings back to nearest edge if non-periodic + inline bool wrap_to_edge(std::vector & ix, std::vector & edge_bin) const + { + bool edge = false; + edge_bin = ix; + for (size_t i = 0; i < nd; i++) { + if (periodic[i]) { + ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined) + edge_bin[i] = ix[i]; + } else if (ix[i] < 0) { + edge = true; + edge_bin[i] = 0; + } else if (ix[i] >= nx[i]) { + edge = true; + edge_bin[i] = nx[i] - 1; + } + } + return edge; + } + + /// \brief Report the bin corresponding to the current value of variable i inline int current_bin_scalar(int const i) const { return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); } + /// \brief Report the flattened bin address corresponding to the current value of all variables + /// and assign first or last bin if out of boundaries + inline int current_bin_flat_bound() const + { + std::vector index = new_index(); + for (size_t i = 0; i < nd; i++) { + index[i] = current_bin_scalar_bound(i); + } + return address(index); + } + /// \brief Report the bin corresponding to the current value of variable i /// and assign first or last bin if out of boundaries inline int current_bin_scalar_bound(int const i) const @@ -447,6 +480,9 @@ template class colvar_grid : public colvarparse { inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const { int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); + + // Wrap bins for periodic dimensions before truncating + if (periodic[i]) bin_index %= nx[i]; if (bin_index < 0) bin_index=0; if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1; return (int) bin_index; @@ -490,6 +526,13 @@ template class colvar_grid : public colvarparse { data[i] = t; } + /// Get the value at the point with linear address i (for speed) + inline T get_value(size_t i) + { + return data[i]; + } + + /// \brief Get the change from this to other_grid /// and store the result in this. /// this_grid := other_grid - this_grid @@ -515,6 +558,7 @@ template class colvar_grid : public colvarparse { has_data = true; } + /// \brief Copy data from another grid of the same type, AND /// identical definition (boundaries, widths) /// Added for shared ABF. @@ -753,6 +797,7 @@ template class colvar_grid : public colvarparse { has_data = true; } + // /// Get the pointer to the binned value indexed by ix // inline T const *value_p (std::vector const &ix) // { @@ -935,10 +980,9 @@ class colvar_grid_count : public colvar_grid } /// \brief Get the binned count indexed by ix from the newly read data - inline size_t const & new_count(std::vector const &ix, - size_t const &imult = 0) + inline size_t const & new_value(std::vector const &ix) { - return new_data[address(ix) + imult]; + return new_data[address(ix)]; } /// Write the current grid parameters to a string @@ -979,7 +1023,7 @@ class colvar_grid_count : public colvar_grid /// Read a grid written by write_multicol(), incrementin if data is true std::istream & read_multicol(std::istream &is, bool add = false); - /// Read a grid written by write_multicol(), incrementin if data is true + /// Read a grid written by write_multicol(), incrementing if add is true int read_multicol(std::string const &filename, std::string description = "grid file", bool add = false); @@ -1017,10 +1061,88 @@ class colvar_grid_count : public colvar_grid has_data = true; } + /// \brief Return the average number of samples in a given "radius" around current bin + /// Really a hypercube of length 2*radius + 1 + inline int local_sample_count(int radius) + { + std::vector ix0 = new_index(); + std::vector ix = new_index(); + + for (size_t i = 0; i < nd; i++) { + ix0[i] = current_bin_scalar_bound(i); + } + if (radius < 1) { + // Simple case: no averaging + if (index_ok(ix0)) + return value(ix0); + else + return 0; + } + size_t count = 0; + size_t nbins = 0; + int i, j, k; + bool edge; + ix = ix0; + // Treat each dimension separately to simplify code + switch (nd) + { + case 1: + for (i = -radius; i <= radius; i++) { + ix[0] = ix0[0] + i; + edge = wrap_detect_edge(ix); + if (!edge) { + nbins++; + count += value(ix); + } + } + break; + case 2: + for (i = -radius; i <= radius; i++) { + ix[0] = ix0[0] + i; + for (j = -radius; j <= radius; j++) { + ix[1] = ix0[1] + j; + edge = wrap_detect_edge(ix); + if (!edge) { + nbins++; + count += value(ix); + } + } + } + break; + case 3: + for (i = -radius; i <= radius; i++) { + ix[0] = ix0[0] + i; + for (j = -radius; j <= radius; j++) { + ix[1] = ix0[1] + j; + for (k = -radius; k <= radius; k++) { + ix[2] = ix0[2] + k; + edge = wrap_detect_edge(ix); + if (!edge) { + nbins++; + count += value(ix); + } + } + } + } + break; + default: + cvm::error("Error: local_sample_count is not implemented for grids of dimension > 3", COLVARS_NOT_IMPLEMENTED); + break; + } + + if (nbins) + // Integer division - an error on the order of 1 doesn't matter + return count / nbins; + else + return 0.0; + } + + /// \brief Return the log-gradient from finite differences /// on the *same* grid for dimension n + /// (colvar_grid_count) inline cvm::real log_gradient_finite_diff(const std::vector &ix0, - int n = 0) + int n = 0, int offset = 0) { cvm::real A0, A1, A2; std::vector ix = ix0; @@ -1028,10 +1150,10 @@ class colvar_grid_count : public colvar_grid // TODO this can be rewritten more concisely with wrap_edge() if (periodic[n]) { ix[n]--; wrap(ix); - A0 = value(ix); + A0 = value(ix) + offset; ix = ix0; ix[n]++; wrap(ix); - A1 = value(ix); + A1 = value(ix) + offset; if (A0 * A1 == 0) { return 0.; // can't handle empty bins } else { @@ -1040,10 +1162,10 @@ class colvar_grid_count : public colvar_grid } } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge ix[n]--; - A0 = value(ix); + A0 = value(ix) + offset; ix = ix0; ix[n]++; - A1 = value(ix); + A1 = value(ix) + offset; if (A0 * A1 == 0) { return 0.; // can't handle empty bins } else { @@ -1054,9 +1176,9 @@ class colvar_grid_count : public colvar_grid // edge: use 2nd order derivative int increment = (ix[n] == 0 ? 1 : -1); // move right from left edge, or the other way around - A0 = value(ix); - ix[n] += increment; A1 = value(ix); - ix[n] += increment; A2 = value(ix); + A0 = value(ix) + offset; + ix[n] += increment; A1 = value(ix) + offset; + ix[n] += increment; A2 = value(ix) + offset; if (A0 * A1 * A2 == 0) { return 0.; // can't handle empty bins } else { @@ -1066,8 +1188,10 @@ class colvar_grid_count : public colvar_grid } } + /// \brief Return the gradient of discrete count from finite differences /// on the *same* grid for dimension n + /// (colvar_grid_count) inline cvm::real gradient_finite_diff(const std::vector &ix0, int n = 0) { @@ -1187,7 +1311,7 @@ class colvar_grid_scalar : public colvar_grid /// Read a grid written by write_multicol(), incrementin if data is true std::istream & read_multicol(std::istream &is, bool add = false); - /// Read a grid written by write_multicol(), incrementin if data is true + /// Read a grid written by write_multicol(), incrementing if add is true int read_multicol(std::string const &filename, std::string description = "grid file", bool add = false); @@ -1259,8 +1383,61 @@ class colvar_grid_scalar : public colvar_grid } } + + /// \brief Return the log-gradient from finite differences + /// on the *same* grid for dimension n + /// (colvar_grid_scalar) + inline cvm::real log_gradient_finite_diff(const std::vector &ix0, + int n = 0, int offset = 0) + { + cvm::real A0, A1, A2; + std::vector ix = ix0; + + // TODO this can be rewritten more concisely with wrap_edge() + if (periodic[n]) { + ix[n]--; wrap(ix); + A0 = value(ix) + offset; + ix = ix0; + ix[n]++; wrap(ix); + A1 = value(ix) + offset; + if (A0 * A1 == 0) { + return 0.; // can't handle empty bins + } else { + return (cvm::logn(A1) - cvm::logn(A0)) + / (widths[n] * 2.); + } + } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge + ix[n]--; + A0 = value(ix) + offset; + ix = ix0; + ix[n]++; + A1 = value(ix) + offset; + if (A0 * A1 == 0) { + return 0.; // can't handle empty bins + } else { + return (cvm::logn(A1) - cvm::logn(A0)) + / (widths[n] * 2.); + } + } else { + // edge: use 2nd order derivative + int increment = (ix[n] == 0 ? 1 : -1); + // move right from left edge, or the other way around + A0 = value(ix) + offset; + ix[n] += increment; A1 = value(ix) + offset; + ix[n] += increment; A2 = value(ix) + offset; + if (A0 * A1 * A2 == 0) { + return 0.; // can't handle empty bins + } else { + return (-1.5 * cvm::logn(A0) + 2. * cvm::logn(A1) + - 0.5 * cvm::logn(A2)) * increment / widths[n]; + } + } + } + + /// \brief Return the gradient of discrete count from finite differences /// on the *same* grid for dimension n + /// (colvar_grid_scalar) inline cvm::real gradient_finite_diff(const std::vector &ix0, int n = 0) { @@ -1288,7 +1465,7 @@ class colvar_grid_scalar : public colvar_grid if (A0 * A1 == 0) { return 0.; // can't handle empty bins } else { - return (A1 - A0) / (widths[n] * 2.); + return cvm::real(A1 - A0) / (widths[n] * 2.); } } else { // edge: use 2nd order derivative @@ -1297,27 +1474,29 @@ class colvar_grid_scalar : public colvar_grid A0 = value(ix); ix[n] += increment; A1 = value(ix); ix[n] += increment; A2 = value(ix); - return (-1.5 * A0 + 2. * A1 - - 0.5 * A2) * increment / widths[n]; + return (-1.5 * cvm::real(A0) + 2. * cvm::real(A1) + - 0.5 * cvm::real(A2)) * increment / widths[n]; } } + /// \brief Return the value of the function at ix divided by its /// number of samples (if the count grid is defined) - virtual cvm::real value_output(std::vector const &ix, - size_t const &imult = 0) const + virtual inline cvm::real value_output(std::vector const &ix, + size_t const &imult = 0) const override { + int s; if (imult > 0) { cvm::error("Error: trying to access a component " "larger than 1 in a scalar data grid.\n"); return 0.; } if (samples) { - return (samples->value(ix) > 0) ? - (data[address(ix)] / cvm::real(samples->value(ix))) : + return ( (s = samples->value(ix)) > 0) ? + (data[address(ix) + imult] / cvm::real(s)) : 0.0; } else { - return data[address(ix)]; + return data[address(ix) + imult]; } } @@ -1325,7 +1504,7 @@ class colvar_grid_scalar : public colvar_grid virtual void value_input(std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0, - bool add = false) + bool add = false) override { if (imult > 0) { cvm::error("Error: trying to access a component " @@ -1334,7 +1513,7 @@ class colvar_grid_scalar : public colvar_grid } if (add) { if (samples) - data[address(ix)] += new_value * samples->new_count(ix); + data[address(ix)] += new_value * samples->new_value(ix); else data[address(ix)] += new_value; } else { @@ -1361,6 +1540,10 @@ class colvar_grid_scalar : public colvar_grid /// \brief Assuming that the map is a normalized probability density, /// calculates the entropy (uses widths if they are defined) cvm::real entropy() const; + + /// \brief Return the RMSD between this grid's data and another one + /// Grids must have the same dimensions. + cvm::real grid_rmsd(colvar_grid_scalar const &other_grid) const; }; @@ -1374,10 +1557,6 @@ class colvar_grid_gradient : public colvar_grid /// should be divided colvar_grid_count *samples; - /// \brief Provide the floating point weights by which each binned value - /// should be divided (alternate to samples, only one should be non-null) - colvar_grid_scalar *weights; - /// Default constructor colvar_grid_gradient(); @@ -1394,6 +1573,13 @@ class colvar_grid_gradient : public colvar_grid /// Constructor from a multicol file colvar_grid_gradient(std::string &filename); + /// Constructor from a vector of colvars and a pointer to the count grid + colvar_grid_gradient(std::vector &colvars, colvar_grid_count *samples_in); + + /// Parameters for smoothing data with low sampling + int full_samples; + int min_samples; + /// Write the current grid parameters to a string std::string get_state_params() const; @@ -1430,25 +1616,25 @@ class colvar_grid_gradient : public colvar_grid cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const; /// Read a grid written by write_multicol(), incrementin if data is true - virtual std::istream & read_multicol(std::istream &is, bool add = false); + std::istream & read_multicol(std::istream &is, bool add = false); /// Read a grid written by write_multicol(), incrementin if data is true - virtual int read_multicol(std::string const &filename, + int read_multicol(std::string const &filename, std::string description = "grid file", bool add = false); /// Write grid in a format which is both human-readable and gnuplot-friendly - virtual std::ostream & write_multicol(std::ostream &os) const; + std::ostream & write_multicol(std::ostream &os) const; /// Write grid in a format which is both human-readable and gnuplot-friendly - virtual int write_multicol(std::string const &filename, + int write_multicol(std::string const &filename, std::string description = "grid file") const; /// Write the grid data without labels, as they are represented in memory - virtual std::ostream & write_opendx(std::ostream &os) const; + std::ostream & write_opendx(std::ostream &os) const; /// Write the grid data without labels, as they are represented in memory - virtual int write_opendx(std::string const &filename, + int write_opendx(std::string const &filename, std::string description = "grid file") const; /// \brief Get a vector with the binned value(s) indexed by ix, normalized if applicable @@ -1474,6 +1660,15 @@ class colvar_grid_gradient : public colvar_grid } } + /// \brief Report the offset between the lower end of the bin and the current value of variable i + inline cvm::real current_offset(int const i) const + { + cvm::real dist = cv[i]->value().real_value - lower_boundaries[i].real_value; + cvm::real offset = std::fmod(dist, widths[i]); + if ( offset < 0. ) offset = widths[i] + offset; + return offset; + } + /// \brief Accumulate the value inline void acc_value(std::vector const &ix, std::vector const &values) { for (size_t imult = 0; imult < mult; imult++) { @@ -1493,28 +1688,89 @@ class colvar_grid_gradient : public colvar_grid samples->incr_count(ix); } - /// \brief Accumulate the gradient based on the force (i.e. sums the - /// opposite of the force) with a non-integer weight - inline void acc_force_weighted(std::vector const &ix, - cvm::real const *forces, - cvm::real weight) { - for (size_t imult = 0; imult < mult; imult++) { - data[address(ix) + imult] -= forces[imult] * weight; - } - weights->acc_value(ix, weight); - } - /// \brief Return the value of the function at ix divided by its /// number of samples (if the count grid is defined) virtual cvm::real value_output(std::vector const &ix, - size_t const &imult = 0) const + size_t const &imult = 0) const override { - if (samples) - return (samples->value(ix) > 0) ? - (data[address(ix) + imult] / cvm::real(samples->value(ix))) : + int s; + if (samples) { + return ( (s = samples->value(ix)) > 0) ? + (data[address(ix) + imult] / cvm::real(s)) : 0.0; - else + } else { return data[address(ix) + imult]; + } + } + + /// Compute the inverse weight corresponding to smoothing factor as in ABF + /// to normalize sums over steps into averages + inline cvm::real smooth_inverse_weight(cvm::real weight) + { + cvm::real fact; + if ( weight <= min_samples ) { + fact = 0.0; + } else if ( weight < full_samples ) { + fact = (weight - min_samples) / (weight * cvm::real(full_samples - min_samples)); + } else { + fact = 1.0 / weight; + } + return fact; + } + + + /// \brief Return the scalar value of the function at ix divided by its + /// number of samples (if the count grid is defined), possibly smoothed + /// by a ramp function going from 0 to 1 between minSamples and fullSamples. + /// Only makes sense if dimension is 1 + virtual inline cvm::real value_output_smoothed(std::vector const &ix, bool smoothed = true) + { + cvm::real weight, fact; + + if (samples) { + weight = cvm::real(samples->value(ix)); + } else { + weight = 1.; + } + + if (smoothed) { + fact = smooth_inverse_weight(weight); + } else { + fact = weight > 0. ? 1. / weight : 0.; + } + + return fact * data[address(ix)]; + } + + /// \brief Obtain the vector value of the function at ix divided by its + /// number of samples (if the count grid is defined), possibly smoothed + /// by a ramp function going from 0 to 1 between minSamples and fullSamples. + inline void vector_value_smoothed(std::vector const &ix, cvm::real *grad, bool smoothed = true) + { + cvm::real weight, fact; + + if (samples) { + weight = cvm::real(samples->value(ix)); + } else { + weight = 1.; + } + + if (smoothed) { + fact = smooth_inverse_weight(weight); + } else { + fact = weight > 0. ? 1. / weight : 0.; + } + + cvm::real *p = &(data[address(ix)]); + + // Appease Clang analyzer, which likes to assume that mult is zero + #ifdef __clang_analyzer__ + assert(mult > 0); + #endif + + for (size_t imult = 0; imult < mult; imult++) { + grad[imult] = fact * p[imult]; + } } /// \brief Get the value from a formatted output and transform it @@ -1523,11 +1779,11 @@ class colvar_grid_gradient : public colvar_grid virtual void value_input(std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0, - bool add = false) + bool add = false) override { if (add) { if (samples) - data[address(ix) + imult] += new_value * samples->new_count(ix); + data[address(ix) + imult] += new_value * samples->new_value(ix); else data[address(ix) + imult] += new_value; } else { @@ -1541,31 +1797,26 @@ class colvar_grid_gradient : public colvar_grid /// Compute and return average value for a 1D gradient grid - inline cvm::real average() + inline cvm::real average(bool smoothed = false) { - size_t n = 0; - if (nd != 1 || nx[0] == 0) { return 0.0; } cvm::real sum = 0.0; - std::vector ix = new_index(); - if (samples) { - for ( ; index_ok(ix); incr(ix)) { - if ( (n = samples->value(ix)) ) - sum += value(ix) / n; - } - } else { - for ( ; index_ok(ix); incr(ix)) { - sum += value(ix); - } + for (std::vector ix = new_index(); index_ok(ix); incr(ix)) { + sum += value_output_smoothed(ix, smoothed); } + return (sum / cvm::real(nx[0])); } + /// \brief Return the RMSD between this grid's data and another one + /// Grids must have the same dimensions. + cvm::real grid_rmsd(colvar_grid_gradient const &other_grid) const; + /// \brief If the grid is 1-dimensional, integrate it and write the - /// integral to a file + /// integral to a file (DEPRECATED by the integrate_potential class) void write_1D_integral(std::ostream &os); }; @@ -1590,12 +1841,16 @@ class integrate_potential : public colvar_grid_scalar integrate_potential(colvar_grid_gradient * gradients); /// \brief Calculate potential from divergence (in 2D); return number of steps - int integrate(const int itmax, const cvm::real & tol, cvm::real & err); + int integrate(const int itmax, const cvm::real & tol, cvm::real & err, bool verbose = true); /// \brief Update matrix containing divergence and boundary conditions /// based on new gradient point value, in neighboring bins void update_div_neighbors(const std::vector &ix); + /// \brief Update matrix containing divergence and boundary conditions + /// called by update_div_neighbors and by colvarbias_abf::adiabatic_reweighting_update_gradient_pmf + void update_div_local(const std::vector &ix); + /// \brief Set matrix containing divergence and boundary conditions /// based on complete gradient grid void set_div(); @@ -1606,6 +1861,63 @@ class integrate_potential : public colvar_grid_scalar add_constant(-1.0 * minimum_value()); } + /// \brief Return the value at the upper edge/corner of the given bin, + /// linearly interpolated (averaged) between 2^nd neighboring bins + /// Used to interpolate to a grid shifted by a half-bin + /// e.g. from PMF to gradient grid + /// should not be used on upper edges of non-periodic grids + /// (that should never be necessary, since those are extra bins with respect + /// to underlying gradient grid) + inline cvm::real value_corner_interp(std::vector const &ix, + size_t const &imult = 0) const + { + cvm::real sum = 0; + std::vector b(ix); + gradients->wrap_to_edge(b, b); + // Fast implementation for small dimensions + if (nd == 1) { + b[0]++; + wrap(b); + return 0.5 * (value(ix, imult) + value(b, imult)); + } + + if (nd == 2) { + sum = value(b, imult); + b[0]++; wrap(b); sum += value(b, imult); + b[1]++; wrap(b); sum += value(b, imult); + b[0]--; wrap(b); sum += value(b, imult); + return 0.25 * sum; + } + + if (nd == 3) { + sum = value(b, imult); + b[0]++; wrap(b); sum += value(b, imult); + b[1]++; wrap(b); sum += value(b, imult); + b[0]--; wrap(b); sum += value(b, imult); + b[2]++; wrap(b); sum += value(b, imult); + b[0]++; wrap(b); sum += value(b, imult); + b[1]--; wrap(b); sum += value(b, imult); + b[0]--; wrap(b); sum += value(b, imult); + return 0.125 * sum; + } + + // Arbitrary dimension + const int max = 1 << nd; + for (int count = 0; count < max; count++) { + count++; + // k-th bit of count is the increment of the k-th coordinate + for (size_t k = 0; k < nd; k++) { + b[k] = ix[k] + (count & 1 << k); + } + wrap_to_edge(b, b); sum += value(b, imult); + } + return sum * (1. / max); + } + + /// \brief Flag requesting the use of a smoothed version of the gradient (default: false) + bool b_smoothed; + + protected: // Reference to gradient grid @@ -1616,10 +1928,6 @@ class integrate_potential : public colvar_grid_scalar // std::vector inv_lap_diag; // Inverse of the diagonal of the Laplacian; for conditioning - /// \brief Update matrix containing divergence and boundary conditions - /// called by update_div_neighbors - void update_div_local(const std::vector &ix); - /// Obtain the gradient vector at given location ix, if available /// or zero if it is on the edge of the gradient grid /// ix gets wrapped in PBC diff --git a/src/colvargrid_def.h b/src/colvargrid_def.h index 92861f43b..fa6531271 100644 --- a/src/colvargrid_def.h +++ b/src/colvargrid_def.h @@ -317,6 +317,9 @@ std::istream & colvar_grid::read_multicol(std::istream &is, bool add) for (size_t i = 0; i < nd; i++ ) { if ( !(is >> x) ) end_of_file = true; bin[i] = value_to_bin_scalar(x, i); + // if x is out of bounds and we are using PBC, wrap it + // Ignore out of bounds points in non-PBC + wrap_detect_edge(bin); } if (end_of_file) break; diff --git a/src/colvarscript_commands_bias.h b/src/colvarscript_commands_bias.h index 420bbabcc..5cecfcaca 100644 --- a/src/colvarscript_commands_bias.h +++ b/src/colvarscript_commands_bias.h @@ -9,7 +9,7 @@ CVSCRIPT(bias_bin, - "Get the current grid bin index (1D ABF only for now)\n" + "Get the current grid bin index (flattened if more than 1d)\n" "bin : integer - Bin index", 0, 0, "", @@ -17,6 +17,8 @@ CVSCRIPT(bias_bin, return COLVARS_OK; ) +/// This is deprecated in the context of mwABF; no other known uses +/// But removing it may break user scripts CVSCRIPT(bias_bincount, "Get the number of samples at the given grid bin (1D ABF only for now)\n" "samples : integer - Number of samples", @@ -36,6 +38,25 @@ CVSCRIPT(bias_bincount, return COLVARS_OK; ) +CVSCRIPT(bias_local_sample_count, + "Get the number of samples around the current bin" + "samples : integer - Number of samples", + 0, 1, + "radius : integer - Sum over radius bins around current bin", + int radius = 0; + char const *arg = + script->obj_to_str(script->get_bias_cmd_arg(0, objc, objv)); + if (arg) { + std::string const param(arg); + if (!(std::istringstream(param) >> radius)) { + script->add_error_msg("local_sample_count: error parsing radius"); + return COLVARSCRIPT_ERROR; + } + } + script->set_result_str(cvm::to_str(this_bias->local_sample_count(radius))); + return COLVARS_OK; + ) + CVSCRIPT(bias_binnum, "Get the total number of grid points of this bias (1D ABF only for now)\n" "Bins : integer - Number of grid points", diff --git a/vmd/cv_dashboard/cv_dashboard.tcl b/vmd/cv_dashboard/cv_dashboard.tcl index 7d64ee82b..6f5dee119 100644 --- a/vmd/cv_dashboard/cv_dashboard.tcl +++ b/vmd/cv_dashboard/cv_dashboard.tcl @@ -120,6 +120,22 @@ Please upgrade to VMD 1.9.4 alpha or later." } } + # Warn about fake PBC dimensions due to non-periodic NAMD: (1, 1, 1) + set abc {0 0 0} + catch { set abc [molinfo $::cv_dashboard::mol get {a b c}]} + if { [vecdist $abc {1 1 1}] < 1e-10 } { + set answer [tk_messageBox -icon warning -type okcancel -title "Warning: unlikely periodic box lengths"\ + -message "The periodic box lengths are (1, 1, 1), which can be the output of a\ + non-periodic simulation in NAMD. Overwrite with (0, 0, 0) to make Colvars detect as non-periodic? + The command line is: + package require pbctools + pbc set {0 0 0} -all -molid $::cv_dashboard::mol"] + if { $answer == "ok" } { + package require pbctools + pbc set {0 0 0} -all -molid $::cv_dashboard::mol + } + } + if {[winfo exists .cv_dashboard_window]} { wm deiconify .cv_dashboard_window ::cv_dashboard::change_track_frame ;# Restart tracking frames when re-opening