pymbar.utils.ParameterError: Warning: Should have \sum_n W_nk = 1

[ThanhGromacs]

Hi everyone,
I just finished 2 sets of Gromacs simulation for my system (RNA+ligand). The first set has only ligand (20 runs). The second set has both RNA and ligand (30 runs). Each simulation is run for 100ns. The force field is AMBER14sb.
The production.mdp for RNA+ligand system:

;====================================================
; Production simulation
;====================================================

;----------------------------------------------------
; RUN CONTROL
;----------------------------------------------------
integrator = sd ; stochastic leap-frog integrator
nsteps = 50000000 ; 2 * 500,000 fs = 1000 ps = 1 ns
dt = 0.002 ; 2 fs
comm-mode = Linear ; remove center of mass translation
nstcomm = 100 ; frequency for center of mass motion removal

;----------------------------------------------------
; OUTPUT CONTROL
;----------------------------------------------------
nstxout = 0 ; don't save coordinates to .trr
nstvout = 0 ; don't save velocities to .trr
nstfout = 0 ; don't save forces to .trr
nstxout-compressed = 1000 ; xtc compressed trajectory output every 1000 steps (2 ps)
compressed-x-precision = 1000 ; precision with which to write to the compressed trajectory file
nstlog = 1000 ; update log file every 2 ps
nstenergy = 1000 ; save energies every 2 ps
nstcalcenergy = 100 ; calculate energies every 100 steps

;----------------------------------------------------
; BONDS
;----------------------------------------------------
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; hydrogens only are constrained
lincs_iter = 1 ; accuracy of LINCS (1 is default)
lincs_order = 4 ; also related to accuracy (4 is default)
lincs-warnangle = 30 ; maximum angle that a bond can rotate before LINCS will complain (30 is default)
continuation = yes ; formerly known as 'unconstrained-start' - useful for exact continuations and reruns

;----------------------------------------------------
; NEIGHBOR SEARCHING
;----------------------------------------------------
cutoff-scheme = Verlet
ns-type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs (default is 10)
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
pbc = xyz ; 3D PBC

;----------------------------------------------------
; ELECTROSTATICS
;----------------------------------------------------
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
ewald_geometry = 3d ; Ewald sum is performed in all three dimensions
pme-order = 6 ; interpolation order for PME (default is 4)
fourierspacing = 0.10 ; grid spacing for FFT
ewald-rtol = 1e-6 ; relative strength of the Ewald-shifted direct potential at rcoulomb

;----------------------------------------------------
; VDW
;----------------------------------------------------
vdw-type = PME
rvdw = 1.0
vdw-modifier = Potential-Shift
ewald-rtol-lj = 1e-3
lj-pme-comb-rule = Geometric
DispCorr = EnerPres

;----------------------------------------------------
; TEMPERATURE & PRESSURE COUPL
;----------------------------------------------------
tc_grps = System
tau_t = 1.0
ref_t = 300
pcoupl = Parrinello-Rahman
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2 ; time constant (ps)
ref_p = 1.0 ; reference pressure (bar)
compressibility = 4.5e-05 ; isothermal compressibility of water (bar^-1)

;----------------------------------------------------
; VELOCITY GENERATION
;----------------------------------------------------
gen_vel = no ; Velocity generation is off

;----------------------------------------------------
; FREE ENERGY CALCULATIONS
;----------------------------------------------------
free-energy = yes
couple-moltype = Protein_chain_B
couple-lambda0 = vdw-q
couple-lambda1 = none
couple-intramol = yes
separate-dhdl-file = yes
sc-alpha = 0.5
sc-power = 1
sc-sigma = 0.3
init-lambda-state = 0
bonded-lambdas = 0.0 0.01 0.025 0.05 0.075 0.1 0.2 0.35 0.5 0.75 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.25 0.5 0.75 1.0 1.00 1.0 1.0 1.0 1.0 1.0 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
vdw-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.05 0.1 0.2 0.3 0.4 0.5 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1.0
nstdhdl = 100
calc-lambda-neighbors = -1

However, using the alchemical_analysis.py to calculate binding free energy on dhdl.xvg files, I encounter an error. (Do I need to modify any file before using alchemical_analysis.py?)

Number of correlated and uncorrelated samples:

State N N_k N/N_k

WARNING: Only 21 uncorrelated samples found at lambda number 0; proceeding with analysis using correlated samples...
0 21 500001 24216.13
1 278 278 1803.88
2 89 89 5656.55
3 1457 1457 343.32
4 97 97 5159.33
5 625 625 800.35
6 2843 2843 175.92
7 3811 3811 131.22
8 562 562 890.48
9 3628 3628 137.84
10 7770 7770 64.36
11 10706 10706 46.70
12 10799 10799 46.30
13 6748 6748 74.10
14 2781 2781 179.81
15 30908 30908 16.18
16 20323 20323 24.60
17 11885 11885 42.07
18 3434 3434 145.62

Estimating the free energy change with MBAR...
Traceback (most recent call last):
File "/Users/thanhle/Desktop/softwares/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1288, in
main()
File "/Users/thanhle/Desktop/softwares/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 1264, in main
Deltaf_ij, dDeltaf_ij = estimatewithMBAR(u_kln, N_k, P.relative_tolerance, regular_estimate=True)
File "/Users/thanhle/Desktop/softwares/alchemical-analysis/alchemical_analysis/alchemical_analysis.py", line 328, in estimatewithMBAR
(Deltaf_ij, dDeltaf_ij, theta_ij ) = MBAR.getFreeEnergyDifferences(uncertainty_method='svd-ew', return_theta = True)
File "/Users/thanhle/anaconda2/lib/python2.7/site-packages/pymbar-3.0.3-py2.7-macosx-10.6-x86_64.egg/pymbar/mbar.py", line 469, in getFreeEnergyDifferences
np.exp(self.Log_W_nk), self.N_k, method=uncertainty_method)
File "/Users/thanhle/anaconda2/lib/python2.7/site-packages/pymbar-3.0.3-py2.7-macosx-10.6-x86_64.egg/pymbar/mbar.py", line 1496, in _computeAsymptoticCovarianceMatrix
check_w_normalized(W, N_k)
File "/Users/thanhle/anaconda2/lib/python2.7/site-packages/pymbar-3.0.3-py2.7-macosx-10.6-x86_64.egg/pymbar/utils.py", line 360, in check_w_normalized
(firstbad, column_sums[firstbad], np.sum(badcolumns)))
pymbar.utils.ParameterError: Warning: Should have \sum_n W_nk = 1. Actual column sum for state 4 was 1.050035. 15 other columns have similar problems

Why are N and N_k columns the same starting from row 2?
I am not quite sure why I got the error. This same error is found in 2 different systems.
Please advice.
Thanks,
Thanh Le

[mrshirts]

The most likely reason for not having the sum equal to zero is that the overlap is poor, and MBAR can't solve it any better.

Try running with alchemical_analysis.py with MBAR deselected in the methods you use. How consistent are TI and BAR? If there are values they don't agree very well, there is likely poor overlap.

[ThanhGromacs]

Thank you for a quick reply,
This might be a foolish question. How does one deselect MBAR? is it this option "free-energy = yes"? If it is correct, I will just rerun all 50 simulations with "free-energy = no"?
Thanks

[mrshirts]

No, this is entirely alchemical_analysis.py. Open the script, and remove MBAR from the list of methods that it loops over. If you don't know what the script is doing at a basic level, you should definitely read it over and figure it out . . .

[ThanhGromacs]

Thank for a quick reply again,
I will do just that.
I really appreciate your help!!!

[ThanhGromacs]

I just finished running the calculation.
Please take a look
TI (kcal/mol)
Coulomb: 351.969 +- 1.623
vdWaals: 40.043 +- 0.150
TOTAL: 392.012 +- 1.630

BAR (kcal/mol)
Coulomb: 388.274 +- 1.400
vdWaals: 41.158 +- 0.151
TOTAL: 429.433 +- 1.408

The difference between 2 methods is about 37 kcal/mol, roughly 9%. Is it too high of a difference?
However TI-CUBIC is very close to BAR.
Coulomb: 381.091 +- 1.366
vdWaals: 39.846 +- 0.148
TOTAL: 420.937 +- 1.374

[mrshirts]

Did you look at the graphical output of alchemical-analysis.py as described by https://doi.org/10.1007/s10822-015-9840-9? That will provide some diagnosis as to where you might need some additional lambda states. just from here, you can see the main issue is lack of agreement in Coulombic energies. 38 or even 7 kcal/mol error is probably not going to provide a good comparison to experiment.

[ThanhGromacs]

Before adding more lambda states to coulomb row, I have a question regarding the result generated by mbar for my complex (RNA+ligand). As you can see from the prod.mdp file, there are 30 simulations for this complex. The weird thing is, I got 0 kcal/mol from states 0 to 10 (I am not sure why they are 0. Comparing these numbers to the wiki tutorial, the trend is similar-not exact) , large numbers from stages 10 to 13 (Coulomb energy decreasing in stages), a large dip for 13-14 (from 75%-100% off), and small number from stages 14-29 (100% off coulomb energy).
Regarding adding more states to coulomb states, what do you think about this new row?
coul-lambdas = 0.0 0.00 0.000 0.00 0.000 0.0 0.00 0.0 0.00 0.0 0.00 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0 1.00 1.0
results.txt
dhdl_TI.pdf
dF_state_long.pdf

Thank you for helping me!