Predicted N-body energies and energy units

[Fe Zn]

Hello,

I am learning about MB-pol and was testing the examples on Github. I was calculating the N-body terms using an example geometry. I am confused what exactly is being computed. From my understanding, the sum of the N-body contributions should be the total energy (in kcal/mol?)?

This is the nrg file:

https://github.com/paesanilab/MBX/blob/master/examples/PEFs/001_mbpol/C%2B%2B/input.nrg

SYSTEM 3H2O
MOLECULE
MONOMER h2o
 O                 -1.58972425    1.04337922   -0.08780840
 H                 -0.63591971    0.97898520    0.00000000
 H                 -1.90066280    1.74501050   -0.66454990
ENDMON
ENDMOL
MOLECULE
MONOMER h2o
 O                  1.64924507    1.08594656    0.00000000
 H                  2.60878026    1.09587704   -0.02817115
 H                  1.33830653    1.78757784    0.57674150
ENDMON
ENDMOL
MOLECULE
MONOMER h2o
 O                 -0.61315209    2.46976336    2.07005086
 H                  0.34684791    2.46976336    2.07005086
 H                 -0.93360667    3.37469919    2.07005086
ENDMON
ENDMOL
ENDSYS

I am using this json configuration file

{"Note" : "This is a  MBX v1.0.0 configuration file","MBX" : {"box" : [] # gas phase ,"twobody_cutoff"   : 500.0 # large cutoff,"threebody_cutoff" : 4.5,"dipole_tolerance" : 1E-08,"dipole_max_it"    : 100,"dipole_method"     : "cg","alpha_ewald_elec" : 0.6,"grid_density_elec" : 2.5,"spline_order_elec" : 6,"alpha_ewald_disp" : 0.6,"grid_density_disp" : 2.5,"spline_order_disp" : 6,"ignore_2b_poly" : [],"ignore_3b_poly" : []} ,"i-pi" : {"localhost" : "localhost3","port" : 34567}}

I run

mb_decomp -v WATER.nrg mbx.json

and get the 1B, 2B, and 3B contributions

0, 2.595209525

1, 1.312228876

2, 1.312228775

0,1, -1.589620258

0,2, 4.165011111

1,2, 2.596461769

0,1,2, 0.3032037847

The sum of which is 10.6947235827.

With the -e flag added,

mb_decomp -v -e WATER.nrg mbx.json

which should print the binding (https://github.com/paesanilab/MBX/blob/master/README.md ) gives

0, 2.595209525

1, 1.312228876

2, 1.312228775

0,1, 2.317818143

0,2, 8.072449411

1,2, 5.22091942

0,1,2, 10.69472358

The 3B term is now 10.69472358, which is exactly the same as the sum from before.

Running the other executable single_point

single_point WATER.nrg mbx.json  also gives Energy= 1.0694723582e+01.

Thank you

[Philip Zhu]

In MBX, the total energy is defined as the binding energy so they should be the same (i.e. the zero of potential in MBX is defined by the configuration in which every monomer is separated infinitely far apart and each monomer is in its optimized geometry). When you run mb_decomp without -e flag, they give you the many-body energy of each sub-cluster. Ex. "0,2, 4.165011111" gives you the 2-body component of monomer[0]&monomer[2]. The sum of everything should be equal to the total binding energy. When you run mb_decomp with the -e flag, the "0,2, 4.165011111" gives you the binding energy of the dimer system monomer[0]&monomer[2]. "0,1,2, 10.69472358" gives you the binding energy of the trimer system monomer[0]&monomer[1]&monomer[2], which should be equal to the total binding energy.

However, I'm not sure I get your question. Your output makes sense to me, and what you said in your question are also true statements that seem to be consistent with your output, so I don't see any issue here. If the above information does not answer your question, can you please clarify?

XZ

[Fe Zn]

Thank you for describing what is being computed, it makes sense now. The part I did not understand before was that the 1-body term is the energy difference between the monomer and an equilibrium geometry.

Just to confirm, the units of the above calculations are in kcal/mol, correct?

Would you be able to provide the total system energy of that equilibrium monomer geometry, as well as what level of theory it was computed at? (If you happen to have the xyz coordinates of the monomer equilibrium geometry that would also be great!). I did not see this discussed in the MB-pol papers I looked at but maybe I missed it.

Thank you very much for the help.

[Philip Zhu]

Yes, unit is in kcal/mol.

The 1-body model is taken from the model developed by Partridge and Schwenke: https://doi.org/10.1063/1.473987. I do not have information about the energy of the 1-body equilibrium geometry, but maybe you can find more information from this paper. 

I do not have the exact equilibrium monomer geometry, but you can optimize an isolated molecule using the optimize executable

[Fe Zn]

Thank you for the link, I went through the paper and the SI. In the paper, they say "r_e and theta_e were fixed at preliminary estimates of the equilibrium geometry.".  I then went into the SI and found the values they used were 0.958649 Angstrom for the OH bond and 104.3475 for the H-O-H angle (in the Fortran code they provided). In that Fortran code, they also provide the energies (in cm-1) they fitted (the c5z array), but the problem is it seems these energies have the equilibrium energy already subtracted because the numbers are not close to what I would expect the energy of a water molecule to be. I could not find that reference energy anywhere in the Fortran code. They do have another file in the SI with the raw data (for example the first data point has an energy of -76.19769220 Ha which is around the energy of a water molecule), but there is no equilibrium geometry reported there either.

I am now a bit confused on how the data used to fit the N-body terms in MB-pol is computed exactly. There must be some reference energy that is subtracted off because without one, there is no way to obtain the 1-body, 2-body, etc. interaction energies that are used to then fit the PIPs used by MB-pol. When collecting the raw data to fit the PIPs, this equilibrium energy must have been used to obtain the interaction energies.

For example for a water dimer geometry, a CCSD(T) calculation would give the total system energy. In MB-pol, the total interaction energy for the dimer is defined to be the sum of two 1-body terms and one 2-body term. To get that 1-body deformation energy, first a CCSD(T) calculation of both of the monomer fragments in that dimer must be computed, after which some equilibrium energy must be subtracted off to get the deformation energy. The 2-body interaction energy is then computed by subtracting the two monomer deformation energies, as well as 2 times the equilibrium energy.

I am attaching an image below with my understanding of how the interaction energies are computed, hopefully I am understanding it correctly. Again, there must be some energy that was subtracted off (E_{ref}) because otherwise there is no way to obtain the interaction energies which are fitted by PIPs.

[Francesco Paesani]

Hello Fe Zn,

You are correct. For example, to calculate the 2-body energy at any level of theory (whether using CCSD(T), DFT, MB-pol, or any force field), one calculates the dimer energy and then subtracts the energy of the two monomers in the same distorted configuration as in the dimer. The same protocol can be applied recursively to any n-body energy, e.g., see Eq. 2 of J. Chem. Theory Comput. 2023, 19, 12, 3551–3566 [https://doi.org/10.1021/acs.jctc.3c00326].

These are two useful references about the MBE. This paper is the first reference to the MBE applied to water: J. Chem. Phys. 53, 4544–4554 (1970) [https://doi.org/10.1063/1.1673986]. This paper provides a nice description of how to use the MBE for calculating high-quality reference energies at the CCSD(T) level: J. Chem. Phys. 135, 224102 (2011) [https://doi.org/10.1063/1.3664730].

We hope this helps, but please don't hesitate to post here if you have additional questions.

All the best,
The MBX Team

[Fe Zn]

Hello,

Thank you for the links, they were very useful. I think I understand what is being computed now for MB-pol, but could you please verify that my understanding is correct. I use the CCSD(T) as the level of theory below but as you said, any level of theory can be applied.

For the 1-body term, the Partridge and Schwenke model is used. That term is the deformation energy with respect to some equilibrium geometry. They refer to these interaction energies as "relaxed" interaction energies in J. Chem. Phys. 135, 224102 (2011) [https://doi.org/10.1063/1.3664730].

For the water dimer, the CCSD(T) total energy of the dimer is calculated. Then, two CCSD(T) calculations are performed to get the total energies of each monomer that make up the dimer. The difference between the total energy of the dimer and the total energies of the monomers is the 2-body interaction energy which is then fitted by PIPs. In total, three CCSD(T) calculations are performed to obtain one 2-body interaction energy.

For the water trimer, the CCSD(T) total energy is calculated, after which three CCSD(T) calculations are done for each of the water monomers in the trimer, as well as three CCSD(T) calculations for each pair of possible dimers in the trimer. The resulting 3-body interaction energy is then fitted by PIPs. In total, 7 CCSD(T) calculations are performed to obtain one 3-body interaction energy.

Just to ask again out of curiosity, have you seen the value of that equilibrium energy mentioned somewhere? I couldn't manage to find a value for it in the Partridge and Schwenke paper. With that equilibrium energy, the total energy of the cluster can be recovered back but the recovered energy might be slightly different than the CCSD(T) total energy that was calculated for the cluster because of the slight differences in level of theory. In any case, that shouldn't really matter because the relative interaction energies are the important part.

[Philip Zhu]

I personally have never looked up for the monomer h2o equilibrium energy used in Partridge and Schwenke's model since this term would be cancelled out in all trainings that we recently perform. However, if you need a number to have an idea, in one of our old scripts, I found a number of -76.37649758 AU for the optimized h2o. I hope this number should be close enough to whatever you might need to play with.

XZ