## 选择性温度积分抽样(SITS)

### 安装配置

将本目录下的所有文件复制粘贴到1.2.6版本的SPONGE目录下，此时会提示是否覆盖（如果没有则弄错了目录），选择覆盖。

然后直接make，在本目录下生成的`SPONGE`文件即为包含SITS的SPONGE模拟文件。

内包含一个测试体系`test_system_ala`，打开后进入其中使用刚生成的SPONGE文件，使用里面对应的mdin文件，正常运行即成功安装SITS模块

### 简介

#### 经典SITS

​		部分稀有事件难以在有限的时间内出现，因此需要使用一些增强抽样算法。温度积分抽样（ITS）是一种以温度进行增强抽样的方式。

​		在NVT系综下，某构象的出现几率服从玻尔兹曼分布：
$$
p(\vec r)\propto{\rm exp}(-\beta_0 U(\vec{r}))\tag{SITS-1}\label{SITS-1}
$$
​		其中$\beta_0=1/{kT}$为系综温度，$U(\vec{r})$为坐标为$\vec{r}$时的能量。

​		通过构建新的有效势能，使得
$$
U_{\rm eff}(\vec r)=-\frac{1}{\beta_0}{\rm ln}\sum_k{n_k{\rm exp}\big(-\beta_kU(\vec r)\big)}\tag{SITS-2}\label{SITS-2}
$$
​		此时得到的分布服从下式
$$
p(\vec r)\propto\sum_k{n_k{\rm exp}\big(-\beta_kU(\vec r)\big)}\tag{SITS-3}\label{SITS-3}
$$
​		其中，$\beta_k$是一系列温度，$n_k$作为指前系数使得每个温度的贡献相等，也即
$$
n_1\sum_{\rm all\ samples}{{\rm exp}\big(-\beta_1U\big)}=n_2\sum_{\rm all\ samples}{{\rm exp}\big(-\beta_2U\big)}\\=\cdots=n_k\sum_{\rm all\ samples}{{\rm exp}\big(-\beta_kU\big)}\tag{SITS-4}\label{SITS-4}
$$


​		动力学过程中的力的变化则为
$$
\vec F_{\rm eff}=-\frac{\part U_{\rm eff}(\vec r)}{\part \vec r}=\frac{\part U_{\rm eff}}{\part U}\vec F=\frac{\sum_k{n_k\beta_k{\rm exp}\big(-\beta_kU\big)}}{\beta_0\sum_k{n_k{\rm exp}\big(-\beta_kU\big)}}\vec F=f_b\vec F\tag{SITS-5}\label{SITS-5}
$$
​		其中，
$$
f_b=\frac{\sum_k{n_k\beta_k{\rm exp}\big(-\beta_kU\big)}}{\beta_0\sum_k{n_k{\rm exp}\big(-\beta_kU\big)}}\tag{SITS-6}\label{SITS-6}
$$
​		最后的统计则需要还原每一份样品的概率，需要乘上$\frac{{\rm exp}(-\beta_0 U)}{\sum_k{n_k{\rm exp}(-\beta_kU)}}$的权重。

​		也可以只对体系的部分能量（$U_1+U_2=U$）进行选择性增强，即为SITS，同样地构建有效势能
$$
U_{\rm eff}(\vec r)=-\frac{1}{\beta_0}{\rm ln}\sum_k{n_k{\rm exp}\big(-\beta_kU_1(\vec r)\big)}+U_2(\vec r)\tag{SITS-7}\label{SITS-7}
$$

$$
p(\vec r)\propto\sum_k{n_k{\rm exp}\big(-\beta_kU_1(\vec r)\big)}{\rm exp}\big(-\beta_0 U_2(\vec r)\big)\tag{SITS-8}\label{SITS-8}
$$

$$
\vec F_{\rm eff}=-\frac{\part U_{\rm eff}(\vec r)}{\part \vec r}=-\bigg(\frac{\part U_{\rm eff}}{\part U_1}\frac{\part U_1}{\part \vec r}+\frac{\part U_{\rm eff}}{\part U_2}\frac{\part U_2}{\part \vec r} \bigg)\\ =\frac{\part U_{\rm eff}}{\part U_1}\vec F_1+\frac{\part U_{\rm eff}}{\part U_2}\vec F_2=\frac{\sum_k{n_k\beta_k{\rm exp}\big(-\beta_kU_1\big)}}{\beta_0\sum_k{n_k{\rm exp}\big(-\beta_kU_1\big)}}\vec F_1+\vec F_2\\=f_b\vec F_1+\vec F_2\tag{SITS-9}\label{SITS-9}
$$

​		还原时同样需要乘上$\frac{{\rm exp}(-\beta_0 U_1)}{\sum_k{n_k{\rm exp}(-\beta_kU_1)}}$的权重。通常能量划分时，将体系划分为AB两部分（AB两部分常常是蛋白质和水，也因此使用字母pw）$U_{\rm pp}$、$U_{\rm ww}$为两部分内部的能量，$U_{\rm pw}$为两部分之间的能量，则可以取$U_1=U_{\rm pp}+f_{\rm pw}U_{\rm pw}$，$U_2=U_{\rm ww}+(1-f_{\rm pw})U_{\rm pw}$。此时
$$
\vec F_{\rm eff}=f_b\vec F_{\rm pp}+(f_bf_{\rm pw}+1-f_{\rm pw})\vec F_{\rm pw}+\vec F_{\rm ww}\tag{SITS-10}\label{SITS-10}
$$


​		在实际操作中，有的时候增强效果不够理想或$n_k$很难收敛，并考虑到能量的平移性和叠加性，可以考虑对有效势能的一些部分进行线性变化，使之变为
$$
U_{\rm eff}(\vec r)=-\frac{1}{\beta_0}{\rm ln}\sum_k{n_k{\rm exp}\big(-\beta_k(aU_1+b)\big)}+U_2+cU_1\tag{SITS-11}\label{SITS-11}
$$

$$
p\propto\sum_k{n_k{\rm exp}\big(-\beta_k (aU_1+b)\big)}{\rm exp}\big(-\beta_0 U_2\big){\rm exp}\big(-\beta_0 cU_1\big)\tag{SITS-12}\label{SITS-12}
$$

$$
\vec F_{\rm eff}=\bigg(\frac{\sum_k{n_k\beta_k{\rm exp}\big(-\beta_k (aU_1+b)\big)}}{\beta_0\sum_k{n_k{\rm exp}\big(-\beta_k(aU_1+b)\big)}}+c\bigg)\vec F_1+\vec F_2=f_b'\vec F_1+\vec F_2\tag{SITS-13}\label{SITS-13}
$$

​		还原时需要乘上$\frac{{\rm exp}(-\beta_0 U_1)}{\sum_k{n_k{\rm exp}(-\beta_k(aU_1+b)}{\rm exp}\big(-\beta_0 cU_1\big)}$的权重。

### IO-经典SITS

#### SITS_nk_traj_file

$n_k$的迭代轨迹文件名。二进制文件，存储$n_k$值，形状为(帧数,`SITS_k_numbers`)。默认值为`SITS_nk_traj_file.dat`。

#### SITS_nk_rest_file

$n_k$的重开文件名。普通文本文件。默认值为`SITS_nk_rest_file.txt`。

#### SITS_nk_in_file

$n_k$的初始文件名。普通文本文件。无默认值。如果该值未设定，$n_k$会被初始化为1.0。如果该值设定，将会读取对应的文件中的值。

#### atom_atribution_in_file

txt 文件。用于指定每个原子属于系统 p/w 的哪一部分，开头一行给出原子数`atom_numbers`，底下`atom_numbers`行每一行给出原子所属的子系统，若为`0`则属于溶质 p，若为`1`则属于溶剂 w。也可以选择下面的 `SITS_protein_numbers` 指定原子所属的子系统。

### 参数-通用

#### SITS_mode

使用的SITS方法。无默认值，必须设定，目前只支持经典 SITS。

- `未设定`：报错。
- `0`：经典 SITS

#### pwwp_enhance_factor

相互作用加强的比例，也即$\eqref{SITS-10}$的$f_{\rm pw}$参数。默认值为`0.5`。

#### SITS_protein_numbers

体系中溶质的原子数，当使用这种方式指定归属时，默认前 `SITS_protein_numbers` 个编号是属于溶质的。也可以使用上面的`atom_atribution_in_file`形式指定原子归属。

### 参数-经典SITS

#### SITS_k_numbers

积分的温度的离散化的份数。默认值为`40`。

#### SITS_T_high

积分的最高温度。

#### SITS_T_low

积分的最高温度。当给定 `SITS_T_high/SITS_T_low` 时温度将依据 `SITS_k_numbers` 在它们之间等间距选取。

#### SITS_T

给出一系列温度，中间用/间隔，例如 `SITS_T = 300/310/320`，这种给定温度的方式和上面用 `SITS_T_high/SITS_T_low` 给定温度的方式选取一个。

#### SITS_pe_a

$\eqref{SITS-11}$中的$a$值。默认值为`1.0`。

#### SITS_pe_b

$\eqref{SITS-11}$中的$b$值。默认值为`0.0`。

#### SITS_fb_shift

$\eqref{SITS-11}$中的$c$值。默认值为`0.0`。

#### SITS_nk_rest

是否从外界读取$n_k$进行初始化，当`SITS_nk_rest = 1`时，必须指定`SITS_nk_rest_file`。

#### SITS_nk_fix

$n_k$是否进行迭代计算。默认值为`0`

- `0`：进行迭代计算。
- 其他值：不进行迭代计算。

#### SITS_record_interval

分子模拟每隔`sits_record_interval`步记录一次能量。默认值为`1`。

#### SITS_update_interval

每记录`sits_update_interval`步能量迭代一次$n_k$。默认值为`100`。

### 如何使用 SITS 加强体系特定部分的能量

SITS 的所有控制信息在代码中存储于 `SITS_INFORMATION` 中，里面留有 `pw_select`接口作为特定能量/力的接口，如果使用者想加强特定部分能量的采样，只需在 `Main_Calculate_Force()` 函数对力/能量进行计算时，将那部分的能量和力依据它们的性质存储到 `pw_select.select_atom_energy[0],pw_select.select_force[0]`（代表 pp 相互作用的能量/力）或者`pw_select.select_atom_energy[1],pw_select.select_force[1]`（代表 pw 相互作用的能量/力）中。例如我想加强溶质/溶质键（bond）的相互作用，只需修改代码：

```c
bond.Bond_Force_With_Atom_Energy_And_Virial(md_info.uint_crd, md_info.pbc.uint_dr_to_dr_cof, md_info.frc, md_info.d_atom_energy, md_info.d_atom_virial);
```

至：

```c
bond.Bond_Force_With_Atom_Energy_And_Virial(md_info.uint_crd, md_info.pbc.uint_dr_to_dr_cof, sits.pw_select.select_force[0], sits.pw_select.select_atom_energy[0], md_info.d_atom_virial);
```

即可，若不想加强这部分的能量，将其该回到上面的形式即可。

**注意**：我们这样做是将体系中的所有 bond 都列入了增强的范畴。如果体系中的水采用的是 TIP3P 力场，并用 Constrain 控制，这么做没有问题；但如果水用的是 flexible TIP3P 力场或其他带成键相互作用的溶剂力场，会使水的成键相互作用也被增强，这不是我们希望看到的。此时需要额外写一个 Bond 的 `Bond_Force_With_Atom_Energy_And_Virial`计算函数 Kernel，利用存储在 `SITS_INFORMATION`中的原子归属信息`atom_sys_mark`进行判断，并将能量和力分别存储到`md_info`和`sits`的相关数组中。

其他部分的能量同理。

