From 70b72d2a6122bde885168b1a2546025efacd8cc8 Mon Sep 17 00:00:00 2001
From: Zhang <bzfzhang@downpour.zib.de>
Date: Tue, 12 Dec 2023 16:43:08 +0100
Subject: [PATCH] Record max norm of applied force

---
 src/colvarbias_abf.cpp | 22 ++++++++++++++++++++++
 src/colvarbias_abf.h   |  3 +++
 src/colvarmodule.cpp   |  3 +++
 src/colvarproxy.cpp    |  2 ++
 4 files changed, 30 insertions(+)

diff --git a/src/colvarbias_abf.cpp b/src/colvarbias_abf.cpp
index 94d41e767..bbbe67768 100644
--- a/src/colvarbias_abf.cpp
+++ b/src/colvarbias_abf.cpp
@@ -29,6 +29,8 @@ colvarbias_abf::colvarbias_abf(char const *key)
     czar_gradients(NULL),
     czar_pmf(NULL),
     last_gradients(NULL),
+    max_force_norm(NULL),
+    max_force_sample(NULL),
     last_samples(NULL)
 {
   colvarproxy *proxy = cvm::main()->proxy;
@@ -190,6 +192,11 @@ int colvarbias_abf::init(std::string const &conf)
   gradients->samples = samples;
   samples->has_parent_data = true;
 
+  max_force_norm = new colvar_grid_scalar(colvars);
+  max_force_sample = new colvar_grid_count(colvars);
+  max_force_norm->samples = max_force_sample;
+  max_force_sample->has_parent_data = true;
+
   // Data for eAB F z-based estimator
   if ( b_extended ) {
     get_keyval(conf, "CZARestimator", b_CZAR_estimator, true);
@@ -295,6 +302,16 @@ colvarbias_abf::~colvarbias_abf()
     gradients = NULL;
   }
 
+  if (max_force_norm) {
+    delete max_force_norm;
+    max_force_norm = NULL;
+  }
+
+  if (max_force_sample) {
+    delete max_force_sample;
+    max_force_sample = NULL;
+  }
+
   if (pmf) {
     delete pmf;
     pmf = NULL;
@@ -356,6 +373,9 @@ int colvarbias_abf::update()
     force_bin = bin;
   }
 
+  if (max_force_sample->index_ok(force_bin)) 
+      max_force_norm->acc_value(force_bin, cvm::proxy->max_atoms_applied_force());
+
   if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
 
     if (update_bias) {
@@ -645,6 +665,8 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool clo
   write_grid_to_file<colvar_grid_count>(samples, prefix + ".count", close);
   write_grid_to_file<colvar_grid_gradient>(gradients, prefix + ".grad", close);
 
+  write_grid_to_file<colvar_grid_scalar>(max_force_norm, prefix + ".mfn", close);
+
   if (b_integrate) {
     // Do numerical integration (to high precision) and output a PMF
     cvm::real err;
diff --git a/src/colvarbias_abf.h b/src/colvarbias_abf.h
index ede514d6e..86323f252 100644
--- a/src/colvarbias_abf.h
+++ b/src/colvarbias_abf.h
@@ -110,6 +110,9 @@ class colvarbias_abf : public colvarbias {
   /// n-dim grid of CZAR pmf (dimension 1 to 3)
   integrate_potential   *czar_pmf;
 
+  colvar_grid_count     *max_force_sample;
+  colvar_grid_scalar *max_force_norm;
+
   inline int update_system_force(size_t i)
   {
     if (colvars[i]->is_enabled(f_cv_subtract_applied_force)) {
diff --git a/src/colvarmodule.cpp b/src/colvarmodule.cpp
index 0b905bd13..3818b18eb 100644
--- a/src/colvarmodule.cpp
+++ b/src/colvarmodule.cpp
@@ -807,6 +807,8 @@ int colvarmodule::calc()
 
   error_code |= analyze();
 
+  proxy->compute_max_atoms_applied_force();
+
   // write trajectory files, if needed
   if (cv_traj_freq && cv_traj_name.size()) {
     error_code |= write_traj_files();
@@ -1924,6 +1926,7 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
        bi++) {
     (*bi)->write_traj(os);
   }
+  os << " " << proxy->max_atoms_applied_force_id() << ' ' << proxy->max_atoms_applied_force();
   os << "\n";
 
   cvm::decrease_depth();
diff --git a/src/colvarproxy.cpp b/src/colvarproxy.cpp
index dbea53265..d288d32c0 100644
--- a/src/colvarproxy.cpp
+++ b/src/colvarproxy.cpp
@@ -7,6 +7,8 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
+#include <iostream>
+
 #include <fstream>
 #include <list>
 #include <utility>