Extreme PMF values for reweighting DBMD

[mallanos]

Hello,
Im trying to use the PyReweighting-2D.py to get the reweighted psi-phi PMF from a DeepBoosted MD simulation. However, the energy values I'm getting are crazy high like 100 Kcal/mol or more.

This is the command Im using to get the input weights file:
awk '!/^#/{if(NR%1==0) print ($8+$7)/(0.001987*300)" " $2 " " ($8+$7)}' ../gamd-0.log > weights.dat

And this is the one I'm using to run the script:
./pyreweighting/reweight-2d.sh 20 10 6 6 coi1_phi-2_psi-2.txt 300

The files gamd-0.log and coi1_phi-2_psi-2.txt are attached

Do you have any idea what's going on?

Thanks

gamd-0.log
coi1_phi-1_psi-1.txt

[MiaoLab20]

Hi Hung, Would you have any thoughts on this? Thanks, Yinglong

…

________________________________ From: Manu Llanos ***@***.***> Sent: Monday, March 18, 2024 3:15:51 PM To: MiaoLab20/pyreweighting ***@***.***> Cc: Subscribed ***@***.***> Subject: [MiaoLab20/pyreweighting] Extreme PMF values for reweighting DBMD (Issue #5) You don't often get email from ***@***.*** Learn why this is important<https://aka.ms/LearnAboutSenderIdentification> Hello, Im trying to use the PyReweighting-2D.py to get the reweighted psi-phi PMF from a DeepBoosted MD simulation. However, the energy values I'm getting are crazy high like 100 Kcal/mol or more. This is the command Im using to get the input weights file: awk '!/^#/{if(NR%1==0) print ($8+$7)/(0.001987*300)" " $2 " " ($8+$7)}' ../gamd-0.log > weights.dat And this is the one I'm using to run the script: ./pyreweighting/reweight-2d.sh 20 10 6 6 coi1_phi-2_psi-2.txt 300 The files gamd-0.log and coi1_phi-2_psi-2.txt are attached Do you have any idea what's going on? Thanks gamd-0.log<https://github.com/MiaoLab20/pyreweighting/files/14641327/gamd-0.log> coi1_phi-1_psi-1.txt<https://github.com/MiaoLab20/pyreweighting/files/14641336/coi1_phi-1_psi-1.txt> — Reply to this email directly, view it on GitHub<#5>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ARJGPFSXCHBJQ7DCWTEZRD3YY44OLAVCNFSM6AAAAABE4GIRHWVHI2DSMVQWIX3LMV43ASLTON2WKOZSGE4TGMBWGQZTKOI>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[hungd238]

To fix your problem, there are two different ways:
(1) You can run more rounds of DBMD equilibration with more steps each round. From what I observed, as the system became more equilibrated, the boost potentials became smaller. DBMD also gave out the averages and standard deviations of boost potentials (along with Vmin and Vmax) after each equilibration round so that users could choose the values they felt reasonable for their simulations.

In the DBMD manuscript, we mentioned that the smaller the nbEndLayers, the larger the boost potentials: https://github.com/MiaoLab20/DBMD/blob/main/deeplearningmodel.py. To maximize the acceleration, for explicit solvent simulations, DBMD automatically uses one nbEndLayer. For implicit solvent simulations, DBMD automatically uses three nbEndLayers to lower the boost potentials. So,
(2) (Not recommended) You can modify lines 34 and 35 of https://github.com/MiaoLab20/DBMD/blob/main/deeplearningmodel.py to use three nbEndLayers (nbEndLayers = 3) to lower the boost potentials.

[mallanos]

Thank you for your quick response!

I plotted the boost potential for the equilibration rounds, but they don't seem to be decreasing over time (Fig attached).
This is the configuration I'm using right now, do you think it looks reasonable for a 4 AA peptide? Should I increase the prep-equi too or just the equilibration. Also, the values chosen for refEP_factor, refED_factor could also be related to the large boot potentials?

Thanks

parmFile = "system.prmtop"
crdFile = "system.rst7"
simType = "explicit" # "explicit", "protein.implicit", "RNA.implicit"

nbCutoff = 10.0 # if "explicit"
temperature = 300

ntcmd = 1000000 # 2ns
cmdRestartFreq = 100

ncycebprepstart, ncycebprepend = 0, 3
ntebpreppercyc = 500000 #1ns
ebprepRestartFreq = 100

ncycebstart, ncycebend = 0, 2
ntebpercyc = 2500000 #5ns
ebRestartFreq = 100

ncycprodstart, ncycprodend = 0, 3
ntprodpercyc = 25000000 #50ns
prodRestartFreq = 500

refEP_factor, refED_factor = 0.05, 0.05 # between 0 and 1

[hungd238]

Thank you for reminding me about refEP_factor and refED_factor. I completely forgot about them. For such small systems, I usually have them at 0. And yes, you may want to run longer pre-equilibration as well to better equilibrate the systems before entering the equilibration stage.

[mallanos]

I've tried to increase the number and lenght of both pre-equilibration and equilibration, but I still get very high boosting potentials. In fact, I noticed that during the pre-equilibration rounds the boost is augmenting each round, so I ended up running only 1 round of pre-equilibration.
This is with nbEndLayers = 3, with nbEndLayers = 1 were even higher.

Then, during the equilibration phase the boosting is still quite high.

The dV vs p(dV) plot seems reasonable, at least to me:

I used the last restart to run production simulations. I simulated 50ns and when I reweight the simulations I still get energies around 20 Kcal/mol

do you have any thougths on what could be happening?

Thanks!

[hungd238]

I have looked through the DBMD codes again, but I found no bug in the codes. It was quite surprising to me that you got quite high boosts during the pre-equilibration stage for the total boost. I can only think of two things left:

(1) Perhaps you may need to run longer cMD before adding some boosts during the pre-equilibration stage. Since the boosts actually reduced the longer you ran the equilibration, I do feel like you may need longer cMD before adding the boosts (so that the system could get to the equilibrium and the gap between Vmin and Vmax would be reduced). And 7.45 on average was actually quite reasonable boost potentials. The highest energy being 20kcal/mol is actually quite low compared to what we usually observe after cumulant expansion to the 2nd order (even for GaMD).

(2) What force field parameter sets do you use for your simulations?

[mallanos]

I could try running longer cMD, currently I'm doing 5ns because the system was run for 250ns before, and I'm using the last snap as input.
I thought 20 Kcal was too much, as most PMFs I see don't go beyond 10kcal-ish. The thing that worries me the most is the fuzziness of the PMF, as even with 150ns I still cant get a more detailed description of the FES. How can I improve that?
I'm using Espaloma parameters for my peptide, instead of a classic protein ff. Do you think that could be a problem? I've tested Espaloma in other contexts and seems to perform well

Thanks!!

[hungd238]

There are two things I would like to clarify here:

(1) We usually use a scale of 0-8 kcal/mol to clearly show the low-energy wells in our publications, but the actual maximum free energy after reweighting can go much higher than that. So as long as you get a reasonable looking FES with clear low-energy wells, it is fine.

(2) Currently, GaMD and DBMD in OpenMM are actually only tested for the AMBER force fields (or AMBER like force fields - D.E. Shaw etc) since to add the boost properly, you need to separate the energy of the whole system into separate energy components. This issue has been brought up here: MiaoLab20/gamd-openmm#22, and that's why you need to specify which type of simulation system (simType) you are running in the https://github.com/MiaoLab20/DBMD/blob/main/simParams.py file so that we can get the correct energy components to add the boost. For AMBER and AMBER-like force fields in explicit solvents, you usually get 6-7 energy components:

FF14SB:
0 HarmonicBondForce
1 HarmonicAngleForce
2 PeriodicTorsionForce
3 NonbondedForce
4 CMMotionRemover
5 MonteCarloBarostat

FF19SB:
0 HarmonicBondForce
1 HarmonicAngleForce
2 PeriodicTorsionForce
3 CMAPTorsionForce
4 NonbondedForce
5 CMMotionRemover
6 MonteCarloBarostat

but for CHARMM (another popular force field), you actually get 10, and that could easily mess up the addition of the boost potentials:

CHARMM36:
0 HarmonicBondForce
1 HarmonicAngleForce
2 HarmonicBondForce (Urey-Bradley)
3 PeriodicTorsionForce
4 CustomTorsionForce (Impropers)
5 CMAPTorsionForce
6 NonbondedForce
7 CustomNonbondedForce (NBFIX)
8 CMMotionRemover
9 MonteCarloBarostat

This is more of a problem with the OpenMM simulation engine rather than the methods. But given the interests I have been receiving for DBMD, I will consider developing it further to resolve this issue in the near future. As of now, I have gone to another institution, so currently need to focus on other projects. However, I will get back to DBMD whenever I have time.

With that said, can you please check in the .mdout files to see how many energy component groups you got and what their names were?

[mallanos]

Hi Hung,
Thanks for your suggestions, and sorry for my late reply, I've been trying to reproduce the Chignolin PMF to test that I'm doing things right.
These are the energy terms of the gamd.log:
ntwx,total_nstep,Total-Energy,Dihedral-Energy,Total-Force-Weight,Dihedral-Force-Weight,Total-Boost,Dihedral-Boost,Total-Harmonic-Force-Constant,Dihedral-Harmonic-Force-Constant,Minimum-Total-Energy,Maximum-Total-Energy,Minimum-Dihedral-Energy,Maximum-Dihedral-Energy,Total-Reference-Energy,Dihedral-Reference-Energy

For Chignolin, I've run 2ns prep, 3x5ns equi, and 2x500ns productions (I'm doing a third one now).
I did the reweighting using cutoff=100 and binsize=1.0 for the RMSD vs Rgyr PMF

dV all: avg = 21.13407400247 std = 5.2708003378739745

when I checked the FES generated by the script, it looks like this:

I wanted to re-plot the data so I grabbed the XVG file corresponding to the cumulant expansion to the 2 order:
fname = f"pmf-2D-c2-{sysname}_rgyr_carmsd.dat-reweight-discx{bin_size}-discy{bin_size}.xvg"
and plotted using a grid interpolation, like this:

# Read data from file
df = pd.read_csv(f'{folder}/reweighting/{fname}', delim_whitespace=True, skiprows=5,
                    names=['RMSD','Rgyr','PMF(kcal/mol)'])

Y = df['Rgyr']
X = df['RMSD']
Z = df['PMF(kcal/mol)']

# Create a grid mesh and interpolate Z
x_min, x_max = X.min(), X.max()
y_min, y_max = Y.min(), Y.max()

xi = np.linspace(x_min, x_max, 1000)
yi = np.linspace(y_min, y_max, 1000)
X_grid, Y_grid = np.meshgrid(xi, yi)
Z_grid = griddata((X, Y), Z, (X_grid, Y_grid), method='cubic', rescale=True)

# Scale Zgrid values
Z_grid = rescale_array(Z_grid, 0, 8)

Z_grid[Z_grid > 6] = 8

levels = [0, 0.3, 0.6, 1.2, 1.8, 3.6, 4.2, 5, 6, 7, 8]

# Plot filled contours
cmap = plt.get_cmap('viridis_r')

CFplot = plt.contourf(X_grid, Y_grid, Z_grid, cmap=cmap, levels=levels, 
antialiased=True, extend='neither')
    
linewidths = 0.2
CLplot = plt.contour(X_grid, Y_grid, Z_grid, levels=levels,
colors='k', linestyles='solid', antialiased=True, linewidths=linewidths)

# Add colorbar and labels
cbar = plt.colorbar(CFplot, label='PMF (kcal/mol)')
# Add the contour line levels to the colorbar
cbar.add_lines(CLplot)
plt.ylabel('Radius of Gyration (Å)');        plt.xlabel('RMSD (Å)')

plt.ylim(4,10);    plt.xlim(0, 8)

plt.title(f'{folder} - Rgyr vs RMSD - bin {bin_size}')

plt.tight_layout()
plt.savefig(f'{folder}/reweighting/plots/PMF_{sysname}_all_rgyr_cbrmsd_{bin_size}.png')
plt.show()

This is how the plot looks like now:

I wanted to make sure this is correct, because different ways of plotting the data make a big difference, but I couldn't find any guidance on best practices for doing it.
For example, I noticed that a low energy well appears in the unfolded zone, which depending on how I set the contour levels it could be very similar in energy to the folded state

Thanks for your help!

[hungd238]

In the Miao Lab, we actually used a software called OriginPro (https://www.originlab.com/) to plot the 2D free energy profiles, with all high free energy values above 10 kcal/mol set to 10 kcal/mol first in the resulting reweighted free energy profile. If the simulations are converged, you usually get some free energy values between 0-6 kcal/mol and other remaining free energy values of some numbers that are usually very high like 20+ (set down to 10 kcal/mol). Other reweighting parameters can be found near the end of this paper: https://pubs.acs.org/doi/full/10.1021/acs.jpclett.3c00926. You can try to follow the protocol to see if you can reproduce Figure 3 in the paper (even without Origin Pro).