NAMD parser does not take temperature as argument

[orbeckst]

The NAMD parser #7 (for FEP) does not take T as an argument and does not return u_nk in proper units of kT.

[orbeckst]

I assume for right now that NAMD writes out energies in kcal/mol, based on NAMD 2.13: NAMD configuration parameters : Output Files

The standard units used by NAMD are Angstroms for length, kcal/mol for energy, Kelvin for temperature, and bar for pressure.

@vtlim @mrshirts can you confirm that the energies in NAMD FEP files are in kcal/mol (in particular our test files in the tyr2ala dataset)?

[vtlim]

The energies in the NAMD FEP output files are in kcal/mol. The output file also includes a column for temperature for the present configuration. Could that lead to a difference causing the assertion error?

[orbeckst]

Thanks for confirming.

Currently, the temperature from the output file is not used, is it?

I don't think that it would be correct to use the instantaneous temperature for kT – the Boltzmann factor comes from a Lagrange multiplier for system at constant T (i.e., constant average temperature) so this T should be the "target temperature" of the thermostat in the simulation. But maybe @mrshirts @davidlmobley want to correct me here.

[mrshirts]

Correct kT is using the target temperature. Instantaneous "temperature" is a quantity that averages to the target temperature over a long enough time (or should if things are working correctly!), but the instantaneous value is not really relevant to anything.

…

On Mon, Apr 15, 2019 at 12:44 PM Oliver Beckstein ***@***.***> wrote: Thanks for confirming. Currently, the temperature from the output file is not used, is it? I don't think that it would be correct to use the instantaneous temperature for kT – the Boltzmann factor comes from a Lagrange multiplier for system at constant T (i.e., constant average temperature) so this T should be the "target temperature" of the thermostat in the simulation. But maybe @mrshirts <https://github.com/mrshirts> @davidlmobley <https://github.com/davidlmobley> want to correct me here. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#75 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABATPVD23G3ABOVVML2HFFTPQTEITANCNFSM4HF3C6EQ> .

[davidlmobley]

Agreed, yes, that we don't want to use the instantaneous temperature. But we do want to allow people to use the correct temperature in their analysis.

[vtlim]

@orbeckst I'm trying to understand where the error is coming from. In the Travis message there are the lines:

>       assert X_delta_f[1] == pytest.approx(delta_f, rel=1e-3)
E       assert 10.116843440154604 == 11.0044412521 +- 1.1e-02

where 10.117 is approximately the input value for the test from 6.031 kcal/mol, so is this message saying that the internal test is returning about 11 kT?

[orbeckst]

EDIT: short answer to your question: yes

@vtlim when you look at

alchemlyb/src/alchemlyb/tests/test_fep_estimators.py

Lines 137 to 149 in f095296

	@pytest.mark.parametrize('X_delta_f', [
	(gmx_benzene_coul_u_nk(), 3.044, 0.01640),
	(gmx_benzene_vdw_u_nk(), -3.033, 0.03438),
	(gmx_expanded_ensemble_case_1(), 75.993, 0.11056),
	(gmx_expanded_ensemble_case_2(), 76.009, 0.11220),
	(gmx_expanded_ensemble_case_3(), 76.219, 0.08886),
	(gmx_water_particle_with_total_energy(), -11.675, 0.065055),
	(gmx_water_particle_with_potential_energy(), -11.724, 0.064964),
	(gmx_water_particle_without_energy(), -11.660, 0.064914),
	(namd_tyr2ala(), 6.031269829/kT_NAMD, 0.069813058/kT_NAMD),
	])
	def test_bar(self, X_delta_f):
	self.compare_delta_f(X_delta_f)

and

alchemlyb/src/alchemlyb/tests/test_fep_estimators.py

Lines 98 to 103 in f095296

	def compare_delta_f(self, X_delta_f):
	est = self.cls().fit(X_delta_f[0])
	delta_f, d_delta_f = self.get_delta_f(est)
	assert X_delta_f[1] == pytest.approx(delta_f, rel=1e-3)
	assert X_delta_f[2] == pytest.approx(d_delta_f, rel=1e-3)

you see that

• X_delta_f is the expected value: 6.031269829 was from the test that you wrote and which passed (and which was in kcal/mol because the NAMD parser used to just work with the raw energies in kcal/mol); I divided it by kT_NAMD to compare it to the updated NAMD parser, which is supposed to return values in kT.
• delta_f is the result of running BAR on the input data

The only real change that I made to namd.extract_unk is

alchemlyb/src/alchemlyb/parsing/namd.py

Line 54 in f095296

win_de_arr = beta * np.asarray(win_de)

If something else is needed or if this is wrong then please feel free to change it.

[vtlim]

I think the difference might be BAR should take in reduced energies but the previous answer I provided to the test had the energies in kcal/mol.

[orbeckst]

In principle, all we are doing here is attaching different units. I haven't looked at the BAR code itself but the equations should not really care if you were to set "1 kT = 1 kcal/mol". From the Amber parser and the MBAR case I know that all that happened was a change of units on the end result.

[vtlim]

Might that affect the weighting factors of the BAR algorithm?

[orbeckst]

Did your run your data through alchemical-analysis and if so, what was the answer and what were the units?

[vtlim]

@orbeckst I can try alchemical-analysis on other NAMD tutorial datasets to see what those answers are. The original code I used also took in data directly from the NAMD fepout file so those units would also have been kcal/mol.

EDIT: I read the earlier message too quickly before responding but just to clarify I have not used alchemical-analysis on my data -- I don't think it supports NAMD. I was using a standalone code modified from a labmate.

[orbeckst]

Thanks.

How should we continue here?

[vtlim]

I'm not sure the best approach to move forward at the moment. I tested a few different versions on the NAMD tutorial files and saw a range of outputs.

source                   tutorial          alc-no reduce    alc-reduce kT (300K)   parseFEP (BAR, T=300K)
----------------------------------------------------------------------------------------------------------
tut02_sphere             -77.6 kcal/mol     -77.1184          -129.4463             -77.4451

tut03_tyr2ala (hydrated) +11.7 kcal/mol       7.1821            12.4539               6.8368
tut03_tyr2ala (isolated)  +4.0 kcal/mol       3.1243             5.2648               3.9506
                          
tut04_crown   (bound)    +61.5 kcal/mol      57.6449            97.1911              57.9465
tut04_crown   (free)     +54.1 kcal/mol      53.4546            89.9684              53.4001

[orbeckst]

Is alc-no reduce the current alchemlyb NAMD parser that does not return reduced energies whereas alc-reduce kT is the one from PR #76? Is parseFEP the other code?

EDIT: If my interpretation is correct and results from PR #76 are inconsistent with both the tutorial value and another code then it suggests that there's something wrong with PR #76. (I'd still liked to know if there's anything in BAR that would be problematic when the units change.)

EDIT 2: tut03_tyr2ala (hydrated) is a strange case: it seems to be the only one for which alc-no reduce disagrees substantially from the tutorial value but agrees with parseFEP. @vtlim do you know what BAR code the tutorial and parseFEP use?

[vtlim]

Is alc-no reduce the current alchemlyb NAMD parser that does not return reduced energies whereas alc-reduce kT is the one from PR #76? Is parseFEP the other code?

Yes.

EDIT 2: tut03_tyr2ala (hydrated) is a strange case: it seems to be the only one for which alc-no reduce disagrees substantially from the tutorial value but agrees with parseFEP. @vtlim do you know what BAR code the tutorial and parseFEP use?

That result was strange to me as well, since the tutorial file says that the data is parsed and analyzed by the ParseFEP in plugin in VMD.

[orbeckst]

For a start someone has to go through the BAR equations and figure out what happens when one changes the units.

[mrshirts]

The BAR and MBAR operate on dimensionless quantities: kT must be in the same units as the energies. Otherwise beta*U becomes something different. . . If BOTH the energies and kT change units, then everything is valid, except of course the free energies will be in those different units.

Currently the software modules assume you are passing in the beta*U quantities (already nondimensionalized). So you multiply the results by whatever kT you divided the energies by before you passed them into BAR.

[orbeckst]

I had a look again and once I converted the new NAMD kT parser values in the table in #75 (comment) to kcal/mol, there's general agreement

                     source  tutorial  alc-no reduce  alc-reduce kT (300K)  parseFEP  alc-reduce kcal/mol (300K)
0              tut02_sphere     -77.6       -77.1184             -129.4463  -77.4451                  -77.170865
1  tut03_tyr2ala (hydrated)      11.7         7.1821               12.4539    6.8368                    7.424532
2  tut03_tyr2ala (isolated)       4.0         3.1243                5.2648    3.9506                    3.138670
3     tut04_crown   (bound)      61.5        57.6449               97.1911   57.9465                   57.941566
4      tut04_crown   (free)      54.1        53.4546               89.9684   53.4001                   53.635672

The exception is still the tutorial value for tut03_tyr2ala (hydrated) of 11.7 kcal/mol but the independent in-house parser + VMD parseFEP also gives 6.8368 kcal/mol so I am inclined to trust the alchemical analysis value.

When I run the test (pytest -k TestBAR) I get for tut03_tyr2ala (hydrated) the alc-reduce kT (300K) value of 11.0044 kT = 6.5604 kcal/mol. Not sure why this is different from @vtlim 's value but it is similar to the parseFEP value. For now, I adopt 11.0044 ± 0.10234 kT as the alchemlyb reference value.

With the new reference value the tests should pass.

---

The table was generated from #75 (comment) with the following code:

import pandas as pd
df =  pd.DataFrame({'source': ['tut02_sphere', 'tut03_tyr2ala (hydrated)', 'tut03_tyr2ala (isolated)', 'tut04_crown   (bound)', 'tut04_crown   (free)'], 'tutorial': [-77.6, +11.7, +4.0, +61.5, +54.1], 'alc-no reduce': [-77.1184, 7.1821 , 3.1243, 57.6449 , 53.4546], 'alc-reduce kT (300K)': [-129.4463, 12.4539, 5.2648, 97.1911, 89.9684], 'parseFEP': [-77.4451, 6.8368, 3.9506, 57.9465, 53.4001]})
kB = 1.9872041e-3
kT = kB * 300
df2 = df.assign(**{"alc-reduce kcal/mol (300K)": lambda x: kT*x["alc-reduce kT (300K)"]})
df2