# Sumhills

**Repository Path**: d2denis/sumhills

## Basic Information

- **Project Name**: Sumhills
- **Description**: No description available
- **Primary Language**: Unknown
- **License**: Not specified
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 0
- **Forks**: 0
- **Created**: 2026-05-20
- **Last Updated**: 2026-06-10

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# Sumhill software(cpp version)
This software is a post MD process for free energy surface calculation. 
The metadynamics log file(myhill.log) and scatter_in contains nearly all the basic information except a few common parameters. 
So it can also double-check the mdout.txt output from SPONGE software.
## Background
Collective variables(CV) $\mathbf{s}(\mathbf{R})$ are low-dimensional variables as the functions of the atomic coordinates of the system, which describe the physical behaviour of slow degree of freedoms(DOFs). For protein-ligand binding system, the CVs can be 3-dimensional coordinate of ligand's center of mass(COM):

$$
  \mathbf{s}=\{c_x,c_y,c_z\}, \quad \mathbf{c}= \frac{\sum_L m_L\cdot \mathbf{R}_L}{\sum_L m_L}
$$
Here, ligand atoms' coordinate is also 3-dimensional: $\mathbf{R}_L=\{L_x,L_y,L_z\}$ atom mass is $m_L$, $L$ means the ligand atoms.

- **Well-tempered Metadynamics** 
Well-tempered distribution(since 2008) is the smoother version of original $p_0$:

$$
p^{WT}(s)=p_0(s)^{\frac{1}{\gamma}} 
$$
where well-tempered factor $\gamma\equiv\frac{\Delta T+T}{T}$, bias Temperature $\Delta T=\frac{\gamma-1}{\beta k_B}>0$ is extra the temperature “boost" via the application of bias potential $V_{bias}$.
The well-tempered biasing potential in Metadynamics:

$$
    V_{WT}(s;t)=\sum_i w(V(s'))g(s(t)-s')
    $$
where the weighting factor is scaled according to the previous step:

$$
    w(V(s'))=w_0 e^{-\frac{V(s')}{k_B\Delta T}}=w_0 e^{-\frac{\beta V(s')}{\gamma-1}}
    $$

### FES Calculation in Metadynamics
In the software, two methods are compared: sumhills and reweight.

- Sumhills method:

$$
  F(s) = \lim_{t\to\infty} C(t) - \frac{\gamma }{\gamma-1}V(s;t)
$$
- Reweight method(optional):

$$
     F(s)=C_0-\frac{1}{\beta}\ln\sum_{t=t_{eq}}^\infty\left[\delta(s'(t)-s)\cdot e^{\beta r_{bias}(s'(t))}\right]
$$

 Note that the discrete delta function $\delta(s'(t)-s)\equiv\begin{cases} 0, & s'(t) \notin \text{grid bin} \\ 1 , & s'(t) \in \text{grid bin} \end{cases}$ is too rough, the reweight result from _histogram grid_ need Gaussian kernal density estimation(KDE) to smooth.


3D Metadynamics's biasing potential $V(\mathbf{s};t)$ is time dependent as the Gaussian hills accumulate on the CV space $\mathbf{s}$ :

$$
    V(\mathbf{s};t)=\sum_{\tau=0}^t  G(\mathbf{s,s}'(\tau))
$$
where the Gaussian hill defined as:

$$
G(\mathbf{s,s}'(\tau))\equiv w_{\tau} \exp\left[-\sum_{d=1,2,3}\frac{\left(s^{(d)}-s'^{(d)}(\tau)\right)^2}{2\sigma_i^2}\right]
$$
Here, 3D Gaussian hill is isotropic because the width for 3-dimensions are the same as $\sigma_i$. The weighting factor $w_\tau$ is also the height of Gaussian hill, slowly decade to zero when well-tempered factor $\gamma$ is defined in the metadynamics software.

Note that sum of hill in MD takes up a lot of time because **every time step** update the bias potential & force!!!
Instead, discretized grid representation of Gaussian hill(see ref1) is conventional, what's more, the number of total points to calculate & save $n_{total}$ = $n_x\cdot n_y\cdot n_z$ is fixed and only update when new hill is add to the bias potential rather than growing hill size.

$$
 G(\mathbf{s}_i,\mathbf{s}'(\tau))\equiv w_{\tau} \exp\left[-\sum_{d=1,2,3}\frac{\left(s_i^{(d)}-s'^{(d)}(\tau)\right)^2}{2\sigma_i^2}\right]
$$
For our software, the grid/scatter points to save the Gaussian hill can be any shape and applied with Great restraining Wall(GRW), the `scatter_file` file contain a mask region should be smaller than cubic grid, which number $n_{mask}< n_{total}$ = $n_x\cdot n_y\cdot n_z$ .
So there are only two basic files to read for the FES calculation:

1. The hill file to discribe $s'(\tau),w_\tau$:

   \# of lines = $n_{hills}= n_{steps}/d\tau$ ($d\tau$ for bias update interval)

   \# of columns = 4 or 5(sink/GRW)

2. The scatter file to discribe grid range, list of $s_i$ and even $\sigma_i$(optional) :
   
   \# of lines = $n_{header}+n_{mask}$

   \# of columns = 3 or 4(catheter)

   The header lines contains 5 lines like this
```
 Remark line, anything
 25 	65.8	0.8
 39 	69.6	0.6
 22 	62.8	0.8
 51	51	51	125
```
which means
```cpp
 Remark line, anything
 x_min   x_max   delta_x
 y_min   y_max   delta_y
 z_min   z_max   delta_z
 nx ny nz nmask
```


Table 1: Each columns' meaning in hill and scatter files.
 
|file   |1st column|2nd column|3rd column|4th column|(5th column)|
|  ---- |  -----   | -------- |  -----   |    --    |   -----    |
|hill   |    X     |   Y      | Z   |  $w_\tau$     |  $p_{max}$ |
|scatter|    X     |   Y      | Z   | $\sigma_{ri}$ |            |

### Mask vs catheter
If the region of interest for scatter points is not regular shape as cubic grid, there are more than two ways to discribe it.
Here, we suppose to use mask for the regions that remain 3D shape such as ball, funnel et. al, on the other hand catheter for 1D curve paths.
The mask version of SinkMeta use isotropic Gaussian hill, which means the Gaussian width $\mathbf{\sigma}$ are all the same for 3-dimensional CV space.
In mask version, the scatter points are equidistant as grid, only the region's edge are clipped. Using SinkMeta, CV will sample within the mask region.
Instead of full space **mask**, lower dimensional subspace such as path is also supported, since the thermodynamic functions such as inner energy, free energy are state-dependent.
The 1D path can accelerate the sampling by frequently transition between the initial(0) and final(last) states.
However, 1D path region ignore the entropic barrier that should be consider in full CV space, which also contribute to free energy profile.
The catheter version of SinkMeta use two kinds of sigmas: the $\sigma_s$ along the path line and $\sigma_{ri}$ perpendicular to line.
The perpendicular $\sigma_{ri}$ is path-dependent, which means that the reaction path now become extensable at the entropy dependent domain, just like the percutaneous coronary intervention restores blood flow to the heart muscle using a catheter, miniature balloon.

![Percutaneous coronary intervention](https://hongkongahc.com/wp-content/uploads/2026/04/imgi_3_500.jpg)

## Code structure

- catheter.hpp: header file
```cpp
...
class Catheter
{
public:
    struct Hill {
    double x, y, z;
    double sigma;
    double weight;
    };

    struct ScatterPoint {
    double x, y, z;
    double rot_x,rot_y,rot_z; ///< rotate vector for catheter
    double sigma_r;
    void Set(const coordinate &coor);
    };
...
    struct CVBorder {
    double xmin,xmax,dx;
    double ymin,ymax,dy;
    double zmin,zmax,dz;
    int nx,ny,nz,ntotal;
    using index = vector<int>;
    void Initial(void);
    coordinate Iterate(const int i,const vector<ScatterPoint>& s_pts);
...
};
```
The class `Catheter` contians three structures: `Hill`,`ScatterPoint` and `CVBorder`.

- catheter.cpp: cpp implement file
```cpp
......
int main(int argc, char *argv[]) {
    ///< Initialization, can read from external inputfile.
    Catheter c;
    string input = "sumhill.txt"; // Default using "sumhill.txt"
    if(argc == 2)
    {
       input = string(argv[1]);
    }
    cout << "The input file is " << input << endl;
    if(!c.Initialize(input))
    {
       return 1;
    }
    if(!c.ReadHill())
    {
        return 2;
    }
    if(!c.ReadScatter())
    {
        return 3;
    }
    if(!c.Reweight())//Sumhill())
    {
        return 4;
    }
    if(c.output == "potential")
    {
        if(!c.WritePotential())
        {
            return 5;
        }
    }
    else if (!c.WriteDireclty())
    {
        return 6;
    }


    return 0;
}
```
- run.sh: execute file.
  Contains shell commands for compling&runing the code.

```bash
#!/bin/bash
g++ -o catheter catheter.cpp -std=c++17 -O2
sleep 1
./catheter reweight.txt
```

## Input
- input file: define overall setting such as Gaussian width $\sigma_s$, which is an important parameter in this software.
Here is an example input file(reweight.txt)
```bash
# input file for reweight
hill_file = myhill.log
scatter_file = mc_mask.txt
history_file = history.log
reweight_file = rw.dat
hist_freq = 1
metaout_file = Meta_Potential.dat
sigma = 0.5
gamma = 50  # well-tempered factor
#catheter = 1
cutoff = 50
output = potential
```
- hill_file: such as myhill.log. Record the 3D real time coordinate $\mathbf{s}'(\tau)$ and weighting factor $w_\tau$
- scatter_file: define the shape of $\sigma_{ri}$ and path scatter point coordinate $\mathbf{s}_i$ of index $i$ 
## Output
- history.log: history dependent bias potential $V(s), r_{bias},RCT$, index $i$ and reweight_factor $\beta r_{bias}$, the former 3 can double-check mdout.txt result.
- Meta_Potential.txt: the final potential $V(s;t)$ at time t=nsteps*dt.
- rw.dat(optional/bug): reweight using kernal density estimation(KDE)
## Reference
(1) Pan, Z.; Li, M.; Chen, D.; Yang, Y. I. A Sinking Approach to Explore Arbitrary Areas in Free Energy Landscapes. JACS Au 2025. DOI: 10.1021/jacsau.5c00460.

## Appendix

### A. C(t)

The convergence of metadynamics depend on the C(t) to be gentle, thus the reweighting need to taken from the gentle C(t) part: the last time after several equilibrium.

Well-tempered case
 
$$
 \begin{align}
    C(t) = \frac{1}{\beta}\left[\ln{Z_0(t)}-\ln{Z_V(t)}\right]\\
       = \frac{1}{\beta}\left[\ln\frac{\int \exp\left(\frac{\gamma\beta}{\gamma-1} V_{meta}(s(t))\right) ds}{\int \exp\left(\frac{\beta}{\gamma-1} V_{meta}(s(t)) \right) ds} \right]
 \end{align}
 $$

### B. KDE
KDE is a **non-parametric** way to estimate the probability density function of a dataset. Unlike GMM, which assumes a mixture of Gaussian distributions, KDE estimates the density by averaging over small, local bumps (kernels) placed on each data point. This makes KDE useful for **smoothly estimating complex distributions** without assuming any particular shape.

### C. Rotate vector for catheter
To be continue ...