MBAR matrix from simulation.nc file

[pga043]

Hi,
How can we create a matrix for MBAR analysis from the simulation.nc file? I can read the .nc file where, the energies are stored as "energies[iteration][replica][state]" but, I am not sure which replica to use for each state and how to obtain the energies corresponding to all the other lambdas for each state (for example using frames from state 1 with all the other lambdas). At the end, I would like to compare these results with results that are obtained from [solvent_dag_results].
Thanks for the help!

Here's what I have tried but, in this trial I mistakenly used energies of replicas as energy from all the other lambdas for a particular state.

import numpy as np
import pymbar
from pymbar import *
import netCDF4 as nc

file =  nc.Dataset('simulation.nc')

nstates = file.variables['energies'].shape[2]
nframes = file.variables['energies'].shape[0]

T = 298.15 # kelvin
beta = 1/(0.001987204259*T)


# Initialize an empty matrix
matrix = np.empty((nstates,0))

'''
units: kT
long_name: energies[iteration][replica][state] is the reduced (unitless) energy of replica 'replica' from iteration 'iteration' evaluated at the thermodynamic state 'state'.
'''
for i in range(nstates):
    energies = file.variables['energies'][0:,0:,i].T
    matrix = np.concatenate((matrix, energies), axis=1)


vector = np.empty(nstates)
vector.fill(nframes)

mbar = MBAR(u_kn=matrix, N_k=vector)

result = mbar.compute_free_energy_differences()

# Free energy and uncertainity in kcal/mol
fe = result['Delta_f'][0,1]/beta
unc = result['dDelta_f'][0,1]/beta

[IAlibay]

@pga043 by matrix - do you mean the N states by N states free energy matrix or the MBAR overlap matrix?

[pga043]

Thanks @IAlibay for the help. Here's what I tried but, I get an error during the analysis.

from openmmtools.multistate import MultiStateReporter
from openfe.protocols.openmm_utils.multistate_analysis import MultistateEquilFEAnalysis
from openff.units import unit

reporter = MultiStateReporter("simulation.nc")
ana = MultistateEquilFEAnalysis(reporter=reporter, sampling_method="repex", result_units=unit.kilocalorie_per_mole)
ana.analyzer._unbiased_decorrelated_N_l

Error:
pymbar.utils.ParameterError: Warning: Should have \sum_n W_nk = 1. Actual column sum for state 0 was 0.164269. 11 other columns have similar problems

Also, I would like to add that when I ran this simulation, I got the final relative free energy difference.
Solvent dG: -28.192480789292194 kilocalorie / mole, err 0.37550400513786314 kilocalorie / mole

[IAlibay]

@pga043 which version of pymbar are you currently using?

[pga043]

These are the versions of pymbar and openfe. I installed openfe via conda but then updated it via git.

pymbar 4.0.2
openfe  0.13.0+265.g844d5de

It's quite possible that the simulation.nc file was generated by a previous version of openfe but, I don't know if the information about the version of openfe is kept in some file or not.

[IAlibay]

@pga043 - unfortunately we've found that in some cases we have the following issues with pymbar 4+, we haven't worked out exactly why (our current guess is that the solver being used within openmmtools isn't the right one).

Our best suggestion for now is to stick to pymbar<4 whilst we try to identify the source of this issue.

[pga043]

@IAlibay thank you for the suggestion and your help. I will try again after downgrading pymbar.