GENESIS

We then use two python scripts to plot the time series of repclia IDs and temperatures. If you are not familiar with python and matplotlib, we provide a simple introduction to data analysis with python in tutorial 4.3.

Content of 02.1_plot_rep_exchange.py:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

parm_rep_data = np.loadtxt("./T-REMD_parmID-repID.log")

N_REP = 16 # total number of replicas
n_row = 4 # number of subplot rows
n_col = 4 # number of subplot columns

fig, axes = plt.subplots(n_row, n_col, figsize=(n_col * 4, n_row * 3), constrained_layout=True, sharex=False, sharey=False)
X = parm_rep_data[:, 0] # common X-axis data: time step
for i_ax in range(N_REP):
    m, n = i_ax // n_col, i_ax % n_col
    Y = parm_rep_data[:, i_ax + 1] # Y-axis data: repID
    axes[m, n].plot(X, Y, ".", c="red")
    axes[m, n].set_xticks([20000000 * j for j in range(6)])
    axes[m, n].set_xticklabels([0, 2, 4, 6, 8, 10], fontsize=12)
    axes[m, n].set_xlim(0, 1e8)
    axes[m, n].set_xlabel(r"MD steps ($\times 10^7$)", fontsize=16)
    axes[m, n].set_yticks([j for j in range(1, 17)])
    axes[m, n].set_yticklabels(["", "", "", "4", "", "", "", "8", "", "", "", "12", "", "", "", "16"], fontsize=12)
    axes[m, n].set_ylim(0, 17)
    axes[m, n].set_ylabel("Replica ID", fontsize=16)
    title_text = "Temperature = {0} K".format(i_ax * 10 + 320)
    axes[m, n].set_title(title_text, fontsize=16)

plt.savefig("Replica_exchange_all.png")

Content of 02.2_plot_temperature_exchange.py:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

rep_parm_data = np.loadtxt("./T-REMD_repID-Temperature.log")

N_REP = 16 # total number of replicas
n_row = 4 # number of subplot rows
n_col = 4 # number of subplot columns

fig, axes = plt.subplots(n_row, n_col, figsize=(n_col * 4, n_row * 3), constrained_layout=True, sharex=False, sharey=False)
X = rep_parm_data[:, 0] # common X-axis data: time step
for i_ax in range(N_REP):
    m, n = i_ax // n_col, i_ax % n_col
    Y = rep_parm_data[:, i_ax + 1] # Y-axis data: repID
    color_ax = [i_ax / N_REP, 1, 1 - i_ax / N_REP]
    axes[m, n].plot(X, Y, "-", c=color_ax)
    axes[m, n].set_xticks([20000000 * j for j in range(6)])
    axes[m, n].set_xticklabels([0, 2, 4, 6, 8, 10], fontsize=12)
    axes[m, n].set_xlim(0, 1e8)
    axes[m, n].set_xlabel(r"MD steps ($\times 10^7$)", fontsize=16)
    axes[m, n].set_yticks([j * 10 + 320 for j in range(16)])
    axes[m, n].set_yticklabels(["320", "", "", "", "360", "", "", "", "400", "", "", "", "440", "", "", ""], fontsize=12)
    axes[m, n].set_ylim(315, 475)
    axes[m, n].set_ylabel("Temperature (K)", fontsize=16)
    title_text = "Replica = {0}".format(i_ax + 1)
    axes[m, n].set_title(title_text, fontsize=16)

plt.savefig("Temperature_exchange_all.png")

Now, let’s execute these scripts:

$ ./01_extract_index_temperature_info.sh
$ ./02.1_plot_rep_exchange.py
$ ./02.2_plot_temperature_exchange.py

These commands will create two text files (T-REMD_parmID-repID.log and T-REMD_repID-Temperature.log) and two figures (Replica_exchange_all.png and Temperature_exchange_all.png).

Please refer to tutorial 11.1 for the explanation of the control file 05_qval_residcg_analysis.inp.

Don’t forget to change the path of GENESIS in line 6 of the bash script 05_calculate_qval.sh:

#!/usr/bin/bash

for irep in {1..16}; do
    cp 05_qval_residcg_analysis.inp tmp.inp
    sed -e "s/RUNID/ID$irep/g" -i tmp.inp
    /home/user/GENESIS/bin/qval_residcg_analysis tmp.inp
    rm -f tmp.inp
done

Now, let’s execute the Q-value calculations:

3.6 MBAR analysis

Now let’s carry out the MBAR analysis:

# change directory
$ cd ../6_MBAR
$ ls
mbar.inp

$ /home/user/GENESIS/bin/mbar_analysis mbar.inp | tee mbar.log

3.7 Calculate PMF of Q-value

Finally, we use the GENESIS analysis tool pmf_analysis to get the PMF of Q-value:

# change directory
$ cd ../7_PMF
$ ls
pmf.inp 01_remove_first_line_of_weight_files.sh 02_plot_PMF_Q.py

We use the script 01_remove_first_line_of_weight_files.sh to remove the first line of the “weight” files, to be consistent with the CV (Q-value) files. The script has the following content:

#!/bin/bash
for j in {1..16}; do
    sed "1d" ../6_MBAR/weight$j.dat > weight$j.dat
done

We then calculate the PMF:

$ /home/user/GENESIS/bin/pmf_analysis pmf.inp | tee pmf.log

This command generates a new file, qval.pmf. We can now plot the PMF with the python script 02_plot_PMF_Q.py. The content of 02_plot_PMF_Q.py is:

#!/usr/bin/env python3

import numpy as np
import matplotlib.pyplot as plt

pmf_data = np.loadtxt("./qval.pmf")

fig, ax = plt.subplots(1, 1, figsize=(6, 4), constrained_layout=False, sharex=False, sharey=False)
ax.plot(pmf_data[:, 0], pmf_data[:, 2], "r-")

ax.set_xticks([0.2 * i for i in range(6)])
ax.set_xticklabels(["0.0", "0.2", "0.4", "0.6", "0.8", "1.0"], fontsize=12)
ax.set_xlim(-0.02, 1.02)
ax.set_xlabel("Q value", fontsize=16)

ax.set_yticks([1 * i for i in range(6)])
ax.set_yticklabels([1 * i for i in range(6)], fontsize=12)
ax.set_ylim(-0.2, 4.2)
ax.set_ylabel("PMF (kcal/mol)", fontsize=16)

plt.savefig("qval_pmf.png", dpi=150)

Now let’s plot the PMF:

$ ./02_plot_PMF_Q.py

This figure shows the PMF of Q value at 400K.