# “To renormalize or not to renormalize ?” in the proton-deuteron scattering
calculations

H. Witała  M. Smoluchowski Institute of Physics, Faculty of Physics, Astronomy
and Applied Computer Science, Jagiellonian University, PL-30059 Kraków, Poland
J. Golak  M. Smoluchowski Institute of Physics, Faculty of Physics, Astronomy
and Applied Computer Science, Jagiellonian University, PL-30059 Kraków, Poland
R. Skibiński  M. Smoluchowski Institute of Physics, Faculty of Physics,
Astronomy and Applied Computer Science, Jagiellonian University, PL-30059
Kraków, Poland

###### Abstract

We discuss two approaches which, by applying the screening method, permit one
to include the long range proton-proton (pp) Coulomb force in proton-deuteron
(pd) momentum-space scattering calculations. In the first one, based on Alt-
Grassberger-Sandhas (AGS) equation, presented in Phys. Rev. C71, 054005 (2005)
and 73, 057001 (2006), one needs to renormalize elastic scattering amplitude
before calculating observables. In the second treatment, proposed by us in
Eur. Phys. Journal A 41, 369 (2009), 41, 385 (2009), and arXiv:2310.03433
[nucl.th], this renormalization is avoided. For the proton induced deuteron
breakup reaction both approaches require renormalization of the corresponding
transition amplitudes. We derive the basic equations underlying both methods
under the assumption that all contributing partial wave states are included
and explain why in our approach renormalization of the elastic scattering
amplitude is superfluous. We show that in order to take into account in the
screening limit all partial waves it is required that four additional terms,
based on the 3-dimensional and partial-wave projected pp Coulomb t-matrices,
identical for both approaches, must appear in transition amplitudes. We
investigate importance of these terms for elastic pd scattering below the
breakup threshold.

The Hamlet-like question in the title arose when two preprints [1] and [2],
both dealing with the problem how to include the long range proton-proton (pp)
Coulomb force in momentum space pd scattering calculations through a screened
Coulomb interaction, were posted. The arguments presented in [2] show that in
the well established approach of Refs. [3, 4] the interplay of the pp Coulomb
potential and the deuteron bound state pole in the neutron-proton t-matrix
makes renormalization of the elastic scattering transition amplitude necessary
prior to calculating observables. Contrary to that, in our approach presented
in [1, 5, 6], one avoids such renormalization. In the following we explain
similarities and differences of both treatments and provide justification why
the renormalization in our method for elastic scattering is unnecessary. We
also discuss a very important problem, indispensable in any treatment of the
long-range Coulomb force: how to take into account, in addition to partial
waves utilised when solving corresponding three-nucleon (3N) scattering
equations, all higher partial wave states.

Let us start with the well established approach of Refs. [3, 4] based on the
AGS equation for the pd transition operator $U$ [7, 8]:

$\displaystyle
U\left|{\Phi}\right\rangle=PG_{0}^{-1}\left|{\Phi}\right\rangle+PtG_{0}U\left|{\Phi}\right\rangle~{},$
(1)

where $P$ is defined in terms of transposition operators,
$P=P_{12}P_{23}+P_{13}P_{23}$, $G_{0}$ is the free 3N propagator and $|\Phi>$
is the initial state composed of a deuteron and a momentum eigenstate of the
proton. The t-matrix $t$ is a solution of the 2-body Lippmann-Schwinger (LS)
equation, with the interaction which contains in case of the pp system in
addition to the nuclear part also the Coulomb pp force (assumed to be screened
and parametrized by some parameter $R$). If the state $U|\Phi>$ is known, the
elastic pd scattering amplitude $<\Phi~{}^{\prime}|U|\Phi>$, with
$|\Phi~{}^{\prime}>$ being the final pd state, can be obtained by quadratures
in the standard manner.

In our approach we use the breakup operator T defined as:

$\displaystyle T=tG_{0}U~{}.$ (2)

It fulfills the 3N Faddeev equation which, when nucleons interact with
pairwise forces only, is given by [9, 8]:

$\displaystyle T|\Phi>=tP|\Phi>+tPG_{0}T|\Phi>~{}.$ (3)

The above form of the Faddeev equation ensures that the T operator reflects
directly the properties of the t-matrix. Here the elastic scattering amplitude
is calculated from solutions of (3) by [8, 9]:

$\displaystyle\left\langle{\Phi^{\prime}}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{\Phi^{\prime}}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{\Phi^{\prime}}\right|PT\left|{\Phi}\right\rangle~{},$
(4)

and the transition amplitude for breakup $<\Phi_{0}|U_{0}|\Phi>$ is expressed
in terms of $T\left|{\Phi}\right\rangle$ by [8, 9]

$\displaystyle\left\langle{\Phi_{0}}\right|U_{0}\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{\Phi_{0}}\right|(1+P)T\left|{\Phi}\right\rangle~{},$
(5)

where $|\Phi_{0}>=|\vec{p}\,\vec{q}\,m_{1}m_{2}m_{3}\nu_{1}\nu_{2}\nu_{3}>$ is
the state of three free outgoing nucleons. In the approach based on the AGS
equation the transition amplitude for breakup is given also by Eq. (5) but
with T replaced by U.

The AGS, (1), as well as the Faddeev, (3), equations are solved in the
momentum-space partial-wave basis $|pq\bar{\alpha}>$:

$\displaystyle|pq\bar{\alpha}>\equiv|pq(ls)j(\lambda\frac{1}{2})I(jI)J(t\frac{1}{2})T>~{},$
(6)

where one can differentiate between the partial wave states $|pq\alpha>$ with
total 2N angular momentum $j$ below some value $j_{max}$: $j\leq j_{max}$, in
which the nuclear, $V_{N}$, as well as the pp screened Coulomb interaction,
$V_{c}^{R}$ (in isospin $t=1$ states only), act, and the states $|pq\beta>$
with $j>j_{max}$, for which only the screened Coulomb force $V_{c}^{R}$ is
present in the pp subsystem. Incorporation of the $|pq\beta>$ states is
indispensable due to the long range nature of the pp Coulomb force and the
necessity to perform finally the screening limit $R\to\infty$. In the
following we derive for both approaches the equations in a subspace restricted
to $|pq\alpha>$ states only, which, however, incorporate all contributions
from the complementary subspace of $|pq\beta>$ states. The states $|pq\alpha>$
and $|pq\beta>$ form together a complete system of states (in the following we
use shorthand notation
$\sum\limits_{\alpha}{\int{p^{2}dpq^{2}dq\left|{pq\alpha}\right\rangle\left\langle{pq\alpha}\right|}}\equiv{\left|{\alpha}\right\rangle\left\langle{\alpha}\right|}$):

$\int{p^{2}dpq^{2}dq}(\sum\limits_{\alpha}{\left|{pq\alpha}\right\rangle}\left\langle{pq\alpha}\right|+\sum\limits_{\beta}{\left|{pq\beta}\right\rangle}\left\langle{pq\beta}\right|)={\left|{\alpha}\right\rangle\left\langle{\alpha}\right|}+{\left|{\beta}\right\rangle\left\langle{\beta}\right|}={\rm
I}\\\ ~{}$ (7)

where I is the identity operator.

Let us start with our approach. Projecting Eq. (3) for $T|\Phi>$ on the
$|pq\alpha>$ and $|pq\beta>$ states one gets the following system of coupled
integral equations [1]:

$\displaystyle\left\langle{pq\alpha}\right|T\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}P\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}{\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|}T\left|{\Phi}\right\rangle$
(8) $\displaystyle+$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}{\left|{\beta^{\prime}}\right\rangle\left\langle{\beta^{\prime}}\right|}T\left|{\Phi}\right\rangle~{},$
(9) $\displaystyle\left\langle{pq\beta}\right|T\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\beta}\right|t_{c}^{R}P\left|{\Phi}\right\rangle+\left\langle{pq\beta}\right|t_{c}^{R}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle~{},$
(10)

where $t_{N+c}^{R}$ and $t_{c}^{R}$ are t-matrices generated by the
interactions $V_{N}+V_{c}^{R}$ and $V_{c}^{R}$, respectively.

Inserting $<pq\beta|T|\Phi>$ from (10) into (9) and using (7) one gets:

$\displaystyle\left\langle{pq\alpha}\right|T\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}P\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}t_{c}^{R3d}P\left|{\Phi}\right\rangle$
(11) $\displaystyle-$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|t_{c}^{R}P\left|{\Phi}\right\rangle$
(12) $\displaystyle+$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle$
(13) $\displaystyle+$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}t_{c}^{R3d}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle$
(14) $\displaystyle-$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|t_{c}^{R}PG_{0}\left|{\alpha^{\prime\prime}}\right\rangle\left\langle{\alpha^{\prime\prime}}\right|T\left|{\Phi}\right\rangle~{}.$
(15)

This is a set of coupled integral equations in the space of the
$\left|{{\alpha}}\right\rangle$ states, which exactly incorporates the
contributions of the pp Coulomb interaction from all partial wave states up to
infinity. It is clear that there is a price to pay for taking into account all
states $|pq\beta>$: the necessity to work with the 3-dimensional Coulomb
t-matrix $t_{c}^{R3d}$, obtained by solving the 3-dimensional LS equation
[10].

Presently it is practically impossible to solve Eq. (15) in its full glory.
The reason are drastic amount of computer resources and of computer time
required to calculate the second and the fifth terms with the 3-dimensional
Coulomb t-matrix. Luckily enough, one can rather easily eliminate them at the
expense of increasing the basis of $\left|{{\alpha}}\right\rangle$ states.
Namely, extending the set $\left|{{\alpha}}\right\rangle$ by adding channels
with higher angular momenta, in which only the pp Coulomb interaction is
present, permits one to completely neglect the four terms in (15) due to their
mutual cancellation: the second with the third and the fifth with the sixth
term. The set (15) is then reduced to:

$\displaystyle\left\langle{pq\alpha}\right|T\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|t_{N+c}^{R}P\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|t_{N+c}^{R}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle~{},$
(16)

which is a basic equation in our approach (in [1] called a simplified one). It
has identical structure as so frequently used 3N Faddeev equation for neutron-
deuteron (nd) scattering [9].

To calculate in our approach the elastic scattering transition amplitude one
needs in (4) the second term
$\left\langle{\vec{p}\vec{q}~{}}\right|T\left|{\Phi}\right\rangle$ composed of
low ($\alpha$) and high ($\beta$) partial wave contributions for $T|\Phi>$.
Using the completeness relation (7) one gets:

$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|T\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|{\alpha^{\prime}}\left\rangle\right\langle{\alpha^{\prime}}\left|T\right|{\Phi}\left\rangle+\left\langle{\vec{p}\vec{q}~{}}\right|t_{c}^{R3d}P\left|{\Phi}\right\rangle-\left\langle{\vec{p}\vec{q}~{}}\right|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|t_{c}^{R}P\left|{\Phi}\right\rangle$
(17) $\displaystyle+$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|t_{c}^{R3d}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle-\left\langle{\vec{p}\vec{q}~{}}\right|{\alpha^{\prime}}\left\rangle\right\langle{\alpha^{\prime}}\left|t_{c}^{R}PG_{0}\right|{\alpha^{\prime\prime}}\left\rangle\right\langle{\alpha^{\prime\prime}}\left|T\right|{\Phi}\left\rangle\right.~{}.$
(18)

To account correctly for contributions from $\left|{{\beta}}\right\rangle$
states again four terms are required, two of which contain 3-dimensional
Coulomb t-matrix. The first one,
$\left\langle{\vec{p}\vec{q}~{}}\right|t_{c}^{R3d}P\left|{\Phi}\right\rangle$,
corresponds to the amplitude of the Rutherford point-deuteron pd scattering
and the second one,
$\left\langle{\vec{p}\vec{q}~{}}\right|t_{c}^{R3d}PG_{0}\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|T\left|{\Phi}\right\rangle$,
is a modification of the first one by nucleon-nucleon (NN) interactions.

Now we derive analogous relations in the approach based on the AGS equation.
Projecting (1) on the $|pq\alpha>$ and $|pq\beta>$ states and using shorthand
notation:

$\sum\limits_{\alpha,\tilde{\alpha}}{\int{p^{2}dpq^{2}dqp^{\prime
2}dp^{\prime}\left|{pq\alpha}\right\rangle
t^{\alpha\tilde{\alpha}}(p,p^{\prime};E-\frac{3}{4m}q^{2})G_{0}\left\langle{pq\tilde{\alpha}}\right|}}\equiv\left|{\alpha}\right\rangle
t^{\alpha}G_{0}\left\langle{\alpha}\right|$

one gets the following system of coupled integral equations:

$\displaystyle\left\langle{pq\alpha}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}{\left\langle{\alpha^{\prime}}\right|}U\left|{\Phi}\right\rangle$
(19) $\displaystyle+$
$\displaystyle\left\langle{pq\alpha}\right|P\left|{\beta^{\prime}}\right\rangle
t_{c}^{R{\beta^{\prime}}}G_{0}\left\langle{\beta^{\prime}}\right|U\left|{\Phi}\right\rangle~{},$
(20) $\displaystyle\left\langle{pq\beta}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\beta}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{pq\beta}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}\left\langle{\alpha^{\prime}}\right|U\left|{\Phi}\right\rangle$
(21) $\displaystyle+$
$\displaystyle\left\langle{pq\beta}\right|P\left|{\beta^{\prime}}\right\rangle
t_{c}^{R{\beta^{\prime}}}\left\langle{\beta^{\prime}}\right|P\left|{\Phi}\right\rangle$
(22) $\displaystyle+$
$\displaystyle\left\langle{pq\beta}\right|P\left|{\beta^{\prime}}\right\rangle
t_{c}^{R{\beta^{\prime}}}G_{0}\left\langle{\beta^{\prime}}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}\left\langle{\alpha^{\prime}}\right|U\left|{\Phi}\right\rangle~{}.$
(23)

Inserting $<pq\beta|U|\Phi>$ from (23) into (20) and using (7) one gets
finally:

$\displaystyle\left\langle{pq\alpha}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}{\left\langle{\alpha^{\prime}}\right|}U\left|{\Phi}\right\rangle$
(24) $\displaystyle-$
$\displaystyle\left\langle{pq\alpha}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{c}^{R{\alpha^{\prime}}}\left\langle{\alpha^{\prime}}\right|P\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|Pt_{c}^{R3d}P\left|{\Phi}\right\rangle$
(25) $\displaystyle-$
$\displaystyle\left\langle{pq\alpha}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{c}^{R{\alpha^{\prime}}}G_{0}{\left\langle{\alpha^{\prime}}\right|}P\left|{\alpha^{\prime\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime\prime}}}G_{0}{\left\langle{\alpha^{\prime\prime}}\right|}U\left|{\Phi}\right\rangle$
(26) $\displaystyle+$
$\displaystyle\left\langle{pq\alpha}\right|Pt_{c}^{R3d}G_{0}P\left|{\alpha^{\prime\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime\prime}}}G_{0}{\left\langle{\alpha^{\prime\prime}}\right|}U\left|{\Phi}\right\rangle~{}.$
(27)

This is a set of coupled integral equations in the space spanned by the
$\left|{{\alpha}}\right\rangle$ states, analogous to (18) in our approach.

Again, extending the set $\left|{{\alpha}}\right\rangle$ by adding a finite
number of channels with higher angular momenta, leads to cancellations between
last four terms and set (27) is reduced to the following basic equation for
approach based on AGS equation [3, 4]:

$\displaystyle\left\langle{pq\alpha}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{pq\alpha}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{pq\alpha}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}\left\langle{\alpha^{\prime}}\right|U\left|{\Phi}\right\rangle~{}.$
(28)

To calculate the elastic scattering transition amplitude
$\left\langle{\Phi^{\prime}}\right|U\left|{\Phi}\right\rangle$ one needs
$\left\langle{\vec{p}\vec{q}~{}}\right|U\left|{\Phi}\right\rangle$ composed of
low ($\alpha$) and high ($\beta$) partial wave contributions for $U|\Phi>$.
Employing the completeness relation (7) and Eq. (27) one gets:

$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|U\left|{\Phi}\right\rangle$
$\displaystyle=$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|PG_{0}^{-1}\left|{\Phi}\right\rangle+\left\langle{\vec{p}\vec{q}~{}}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime}}}G_{0}{\left\langle{\alpha^{\prime}}\right|}U\left|{\Phi}\right\rangle$
(29) $\displaystyle-$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{c}^{R{\alpha^{\prime}}}{\left\langle{\alpha^{\prime}}\right|}P\left|{\Phi}\right\rangle+\left\langle{\vec{p}\vec{q}~{}}\right|Pt_{c}^{R3d}P\left|{\Phi}\right\rangle$
(30) $\displaystyle-$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|P\left|{\alpha^{\prime}}\right\rangle
t_{c}^{R{\alpha^{\prime}}}G_{0}{\left\langle{\alpha^{\prime}}\right|}P\left|{\alpha^{\prime\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime\prime}}}G_{0}{\left\langle{\alpha^{\prime\prime}}\right|}U\left|{\Phi}\right\rangle$
(31) $\displaystyle+$
$\displaystyle\left\langle{\vec{p}\vec{q}~{}}\right|Pt_{c}^{R3d}G_{0}P\left|{\alpha^{\prime\prime}}\right\rangle
t_{N+c}^{R{\alpha^{\prime\prime}}}G_{0}{\left\langle{\alpha^{\prime\prime}}\right|}U\left|{\Phi}\right\rangle~{}.$
(32)

Using relation (2) between $U$ and $T$ one finds that indeed amplitudes and
thus also observables are the same in both treatments.

It should be emphasized that only by extending the set of
$\left|{{\alpha}}\right\rangle$ states is it possible to neglect in (15) and
(27) the terms which contain the 3-dimensional Coulomb t-matrices, and to
reduce the problem in both approaches to numerically well treatable equations
(16) and (28). The indication that cancellations takes place is given by
convergence of predictions with respect to the total angular momentum in the
two-nucleon (2N) subsystem $j_{max}$, which defines the set of
$\left|{{\alpha}}\right\rangle$ states. It will be denoted in the following by
$j_{s}\,j_{max}$ with $j_{s}$ being the largest angular momentum in which the
2N interaction acts [1].

It is evident that a correct treatment of the Coulomb force in both approaches
requires inclusion of four additional terms in the elastic (and also breakup)
transition amplitudes (the last four terms in (18) and (32)).

It was shown in [2] (see also [3, 4] and references therein) that in the
treatment based on AGS equation (28) the elastic scattering transition
amplitude acquires in the screening limit $R\to\infty$ an infinitely
oscillating phase factor and must be renormalized before calculating
observables. As a consequence, each term in (32) containing
$U\left|{\Phi}\right\rangle$ has to be renormalized. In our approach we solve
instead of AGS the 3N Faddeev equation (16) for the
$\left\langle{pq\alpha}\right|PT\left|{\Phi}\right\rangle$ states, from which
later elastic scattering transition amplitude is calculated. In this way we
avoid the main source of the oscillating phase factor described in [2] and the
necessity of renormalization of the elastic scattering amplitude.
Additionally, the structure of 3N Faddeev equation guarantees that their
solutions inherit properties from the two-nucleon t-matrices providing thus an
additional argument that renormalization is redundant. Namely, the properties
of t-matrices generated by the screened Coulomb force alone (in the case of
partial wave decomposed t-matrices also those generated by a combination of
Coulomb and nuclear parts) as well as their screening limits were studied
theoretically in the past in numerous papers [11, 12, 13, 14, 15, 16, 17, 18,
19] and later some of these properties were confirmed numerically in [10]. The
most important finding was that such off-shell t-matrices have a well defined
screening limit while the half- and on-shell ones acquire in this limit an
infinitely oscillating phase factor. At the same time, the elastic pd
scattering amplitude gets contributions of
$\left\langle{pq\alpha}\right|T\left|{\Phi}\right\rangle$ states only from the
off-shell region of the Jacobi momenta magnitudes $q$ and $p$ in ($q-p$)
plane:
$\frac{p^{2}}{m}+\frac{3}{4m}q^{2}\neq\frac{3}{4m}q_{max}^{2}=\frac{3}{4m}q_{0}^{2}+E_{d}$,
where $m$ is the nucleon mass, $E_{d}$ is the (negative) deuteron binding
energy, and $q_{0}$ is the magnitude of the relative pd momentum. That off-
shell region of $q-p$ values does not overlap with the ellipse from which
half-on-shell contributions to the breakup reaction come. In Fig. 1 we
exemplify that off-shell part and the separation of the breakup and elastic
scattering regions in the $(q-p)$ plane for the energy of a pd system $E=3.5$
MeV, which is slightly above the breakup threshold and for which both
reactions are possible, and at $E=3.0$ MeV, which is below the breakup
threshold and for which only elastic scattering is allowed. The fact that
elastic pd scattering requires only off-shell solutions of the Faddeev
equations and that the off-shell two-nucleon t-matrices have a well defined
screening limit is the reason why in our method no renormalization of elastic
scattering amplitudes is needed. Contrary to that, the breakup amplitudes
acquire the oscillating phase factor originating from half-shell t-matrices.

In order to compare results of two approaches and check that indeed our method
does not need the renormalization, we applied our approach at a low proton
energy below the breakup threshold, where effects of the pp Coulomb force as
well as contributions of different terms to the elastic scattering amplitude
are expected to be dominant and where also results of the AGS approach are
available at $E=3.0$ MeV [3]. In Fig. 2 we show our predictions obtained with
the AV18 NN potential [20] and $j_{s}3j7$ $\left|{\alpha}\right\rangle$ basis
at $3.0$ MeV compared to existing elastic scattering data for the cross
section and analyzing powers. The red short dashed line show results obtained
with only the first three terms in elastic scattering transition amplitude
(18), which is the approximation used also in Ref. [3]. The red solid lines
are predictions for neutron-deuteron scattering. It is clear that in this
region of energies the Coulomb force effects indeed are large and dominant at
all angles as evidenced by comparing the red solid and short dashed lines. It
is astonishing how good the overall description of tensor analyzing power data
is in spite of their small magnitudes of $\approx 1\%$. The vector analyzing
powers $A_{y}$ and $iT_{11}$ are underestimated by theory what is very well
known in the literature under the name “low energy analyzing power puzzle”.
Even more interesting is the good agreement for practically all shown
observables, with the exception of $A_{y}$ and $iT_{11}$, between our $3.0$
MeV results and the predictions based on the AGS approach, as far as it can be
judged from Fig. 9 of Ref. [3]. This good agreement strongly supports the
statement that both approaches have to provide the same predictions for all
observables and that in our approach renormalization of the elastic scattering
amplitude is indeed superfluous. The differences for $A_{y}$ and $iT_{11}$ can
be very probably traced back to the well known large sensitivity of these
observables to the ${}^{3}P_{j}$ components of the NN interaction [9] and
different dynamics used by us and in [3].

In Fig. 2 we show also by dotted blue lines results with the last term in (18)
included. It is evident that the term
$-\left\langle{\vec{p}\vec{q}~{}}\right|\sum\limits_{\alpha^{\prime}}{\int{\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|}}t_{c}^{R}PG_{0}S\sum\limits_{\alpha^{\prime\prime}}{\int{\left|{\alpha^{\prime\prime}}\right\rangle\left\langle{\alpha^{\prime\prime}}\right|}}T\left|{\Phi}\right\rangle$
is significant at low energies and that it deteriorates good description of
data obtained with the first three terms. In [1] it was shown that at energies
above $\approx 10$ MeV the contribution of that term to elastic scattering
observables is negligible and at $10$ MeV it starts to influence some spin
observables. It is thus unavoidable below the breakup threshold to investigate
how significant are effects of inclusion of the fourth term
$\left\langle{\vec{p}\vec{q}~{}}\right|t_{c}^{R3d}PG_{0}\sum\limits_{\alpha^{\prime}}{\int{\left|{\alpha^{\prime}}\right\rangle\left\langle{\alpha^{\prime}}\right|}}T\left|{\Phi}\right\rangle$
in the elastic scattering transition amplitude. Since the fifth term has a
negative sign and contains partial wave contributions to the Coulomb t-matrix
whose full 3-dimensional form is contained in the fourth term, one would
expect that they would at least partially cancel each other and the inclusion
of the fourth term should restore at least partly the good description of
data.

The computation of the fourth term with the 3-dimensional Coulomb t-matrix
$t_{c}^{R3d}$, can be done according to expressions (D.9), (D.6), and (D.8) of
Ref. [5]. It requires integrations over components of two vectors: over vector
$\vec{q}$ in (D.9), and over $\vec{p}~{}^{\prime}$ or $\vec{q}_{4}$ in (D.6)
or (D.8), respectively. Below the breakup threshold only channels
$\alpha\neq\alpha_{d}$ contribute to (D.6). Since below the breakup threshold
the decomposition (D.7) is superfluous, (D.8) provides the full contribution
from $\alpha_{d}$ channels, obtained by replacing the second part of splitting
(D.7) with the left side of (D.7). The contributions from (D.6) and (D.8) must
be determined numerically and this is the most time consuming part of the
calculations.

In Fig. 2 the indigo crosses show the results obtained with all the terms in
(18) included. As expected the fourth and fifth terms cancel each other to a
large extent and a good description of data for the cross section and tensor
analyzing powers is essentially regained.

To get an idea about the magnitude of the Coulomb force effects for other
elastic scattering observables we show in Figs. 3-6 analogous predictions as
in Fig. 2 but for selected spin correlations (Fig. 3), proton to proton (Fig.
4), proton to deuteron (Fig. 5), and deuteron to proton (Fig. 6) spin transfer
coefficients. The figures reveal a wide spectrum of importance and magnitude
of the Coulomb force effects, dependent on the observable. For most of
observables the effects are large in a wide range of angles, for example for
spin correlations from Fig. 3 and some of spin transfers ($K_{y}^{x}(N-N)$,
$K_{y}^{xz}(N-D)$, $K_{zz}^{y}(D-N)$). For some large effects are restricted
to forward region of angles below $\approx 90^{o}$ ($K_{x}^{x}(N-N)$,
$K_{z}^{x}(N-N)$, $K_{z}^{z}(N-N)$, $K_{x}^{x}(N-D)$, $K_{z}^{x}(N-D)$,
$K_{z}^{z}(N-D)$, $K_{x}^{z}(D-N)$, $K_{z}^{z}(D-N)$). There are some
interesting cases of observables which for the neutron-deuteron scattering
vanish and become nonzero for the proton-deuteron interaction, as for example
the nucleon to nucleon spin transfer coefficient $K_{y}^{x}(N-N)$ shown in
Fig. 4. These nonzero values are due to a large charge independence breaking
of pp and neutron-proton (np) interactions in isospin $t=1$ states, caused by
the Coulomb pp force. In our calculations we used the charge dependent AV18
potentials, taking np and pp NN interactions of this model for the pd and nd
systems. In all isospin $t=1$ states both total isospins of the 3N system
$T=\frac{1}{2}$ and $T=\frac{3}{2}$ were taken into account. Vanishing of the
$K_{y}^{x}(n-n)$ for nd scattering shows that the difference between np and pp
NN AV18 potentials is too weak to induce nonzero values for this observable.

The very interesting and most important effect seen in all figures is that
practically in all cases (large) effects caused by adding the fifth term to
the elastic scattering transition amplitude are removed when including
simultaneously the fourth term. In consequence, it is needless to account for
these terms in elastic scattering amplitude what drastically simplifies and
accelerates determination of the Coulomb force effects.

Summarizing, we have shown that the two discussed approaches which enable to
include the long range Coulomb force in momentum-space pd scattering
calculations by applying a screening method have to provide the same results
for all observables. In each method the cancellation between terms containing
3-dimensional and partial wave decomposed Coulomb t-matrices is decisive for
establishing workable equations, whose structure is identical to the commonly
used equations for neutron-deuteron scattering. Solutions of these equations
together with four additional terms, two of which contain the 3-dimensional
Coulomb t-matrices, permit one to get the elastic scattering (and breakup)
transition amplitudes. In the approach based on the AGS equation it is
unavoidable to perform renormalization of the elastic scattering amplitudes
before calculating observables. In the approach based on the Faddeev equation
such renormalization can be completely avoided. We have shown numerically that
the cancellation of last two terms in elastic scattering transition amplitude
enables one to determine the pp Coulomb force effects in the pd scattering
nearly as easily as to compute observables in neutron-deuteron scattering.

###### Acknowledgements.

This research was supported in part by the Excellence Initiative – Research
University Program at the Jagiellonian University in Kraków. The numerical
calculations were partly performed on the supercomputers of the JSC, Jülich,
Germany.

## References

  * [1] H. Witała, J. Golak, and R. Skibiński, arXiv:2310.03433 [nucl.th].
  * [2] A. Deltuva, arXiv:2311.14605v1 [nucl.th].
  * [3] A. Deltuva, A. C. Fonseca, and P. U. Sauer, Phys. Rev. C71, 054005 (2005).
  * [4] A. Deltuva, A. C. Fonseca, and P. U. Sauer, Phys. Rev. C72, 054004 (2005).
  * [5] H. Witała, R. Skibiński, J. Golak, W. Glöckle, Eur. Phys. Journal A41, 369 (2009).
  * [6] H. Witała, R. Skibiński, J. Golak, W. Glöckle, Eur. Phys. Journal A41, 385 (2009).
  * [7] E. O. Alt, P. Grassberger, W. Sandhas, Nucl. Phys. B2, 167 (1967).
  * [8] W. Glöckle, The Quantum Mechanical Few-Body Problem, Springer Verlag 1983.
  * [9] W. Glöckle, H. Witała, D. Hüber, H. Kamada, J. Golak, Phys. Rep. 274, 107 (1996).
  * [10] R. Skibiński, J. Golak, H. Witała, and W.Glöckle, Eur. Phys. Journal A40, 215 (2009).
  * [11] E. O. Alt, W. Sandhas, and H. Ziegelmann, Phys. Rev. C 17, 1981 (1978).
  * [12] J.C.Y. Chen and A.C. Chen, in Advances of Atomic and Molecular Physics, edited by D. R. Bates and J. Estermann ( Academic, New York, 1972), Vol. 8.
  * [13] L. P. Kok, H. van Haeringen, P hys. Rev. C21, 512 (1980).
  * [14] W.F. Ford, Phys. Rev. 133, B1616 (1964).
  * [15] W.F. Ford, J. Math. Phys. 7, 626 (1966).
  * [16] J.R. Taylor, Nuovo Cimento B23, 313 (1974).
  * [17] M.D. Semon and J.R. Taylor, Nuovo Cimento A26, 48 (1975).
  * [18] L. P. Kok, H. van Haeringen, Phys. Rev. Lett. 46, 1257 (1981).
  * [19] H. van Haeringen, Charged Particle Interactions, Theory and Formulas, (Coulomb Press, Leyden, 1985).
  * [20] R. B. Wiringa, V. G. J. Stoks, R. Schiavilla, Phys. Rev. C51, 38 (1995).
  * [21] S. Shimizu et al., Phys. Rev. C 52, 1193 (1995).

Figure 1: (color online) Regions of the Jacobi momenta $q$ and $p$ values in
$(q-p)$ plane which contribute to the breakup reaction ((red) solid line at
$E=3.5$ MeV, showing ellipse
$\frac{p^{2}}{m}+\frac{3}{4m}q^{2}=\frac{3}{4m}q_{max}^{2}=\frac{3}{4m}q_{0}^{2}+E_{d}$
) and elastic scattering ($<\Phi^{\prime}|PT|\Phi>$ term) (gray highlighted
region) at the incoming nucleon laboratory energy $E=3.0$ and $3.5$ MeV.
Figure 2: (color online) Comparison of data and predictions for the pd
scattering cross section $\frac{d\sigma}{d\Omega}$, proton vector $A_{y}$,
deuteron vector $iT_{11}$ and deuteron tensor $T_{20}$, $T_{21}$, $T_{22}$
analyzing powers. They are shown as functions of a c.m. proton scattering
angle $\Theta_{cm}$ and were calculated at the incoming proton laboratory
energy $E=3.0$ MeV with the approach based on Faddeev equation (16) and
transition amplitude (18). The exponentialy screened Coulomb force ($R=40$ fm,
$n=4$) and the AV18 potential [20] restricted to the $j\leq 3$ partial waves
have been applied. To solve Faddeev equation the set $j_{s}3j7$ of
$\left|{\alpha}\right\rangle$ states was used. The red short dashed lines show
the results when only the first three terms in (18) are taken into account.
The blue dotted lines are predictions when also the fifth term in (18)
($-\left\langle{\vec{p}\vec{q}~{}}\right|{\alpha^{\prime}}\left\rangle\right\langle{\alpha^{\prime}}\left|t_{c}^{R}PG_{0}\right|{\alpha^{\prime\prime}}\left\rangle\right\langle{\alpha^{\prime\prime}}\left|T\right|{\Phi}\left\rangle\right.$)
is included. The pure Coulomb term
$\left\langle{\Phi^{\prime}}|Pt_{c}P|{\Phi}\right\rangle$ was determined using
the screening limit expresion for the off-shell 3-dimensional Coulomb t-matrix
(Eq. (19) in Ref. [1]). The indigo crosses show the results with all terms in
(18) included. The red solid lines are predictions for nd elastic scattering
and green circles represent the pd data from Ref. [21].  Figure 3: (color
online) The same as in Fig. 2 but for selected spin correlation coefficients.
For description of lines see Fig. 2.  Figure 4: (color online) The same as in
Fig. 2 but for selected proton to proton spin transfer coefficients. For
description of lines see Fig. 2.  Figure 5: (color online) The same as in Fig.
2 but for selected proton to deuteron spin transfer coefficients. For
description of lines see Fig. 2.  Figure 6: (color online) The same as in Fig.
2 but for selected deuteron to proton spin transfer coefficients. For
description of lines see Fig. 2.


# MLIMC: Machine learning-based implicit-solvent Monte Carlo

Jiahui Chen1, Weihua Geng2iii Corresponding author, Email<EMAIL_ADDRESS>, and
Guo-Wei Wei1,3,iiiiii Corresponding author, Email<EMAIL_ADDRESS>
1 Department of Mathematics, Michigan State University, MI 48824, USA.
2 Department of Mathematics, Southern Methodist University, Dallas, TX 75275,
USA
3 Department of Biochemistry and Molecular Biology, Michigan State University,
MI 48824, USA.

This work is dedicated to Professor John ZH Zhang on the occasion of his 60th
birthday

###### Abstract

Monte Carlo (MC) methods are important computational tools for molecular
structure optimizations and predictions. When solvent effects are explicitly
considered, MC methods become very expensive due to the large degree of
freedom associated with the water molecules and mobile ions. Alternatively
implicit-solvent MC can largely reduce the computational cost by applying a
mean field approximation to solvent effects and meanwhile maintains the atomic
detail of the target molecule. The two most popular implicit-solvent models
are the Poisson-Boltzmann (PB) model and the Generalized Born (GB) model in a
way such that the GB model is an approximation to the PB model but is much
faster in simulation time. In this work, we develop a machine learning-based
implicit-solvent Monte Carlo (MLIMC) method by combining the advantages of
both implicit solvent models in accuracy and efficiency. Specifically, the
MLIMC method uses a fast and accurate PB-based machine learning (PBML) scheme
to compute the electrostatic solvation free energy at each step. We validate
our MLIMC method by using a benzene-water system and a protein-water system.
We show that the proposed MLIMC method has great advantages in speed and
accuracy for molecular structure optimization and prediction.

Key words: Machine learning, Implicit-solvent Monte Carlo simulation, Poisson-
Boltzmann equation, electrostatics.

##  I Introduction

The determination of protein structures is of paramount importance for
structural biology and macromolecular study. However, not all protein
structures can be determined with available experimental techniques due to
various limitations. Computational methods offer important alternative
approaches for structural determination and optimization [63]. Indeed,
molecular force field models and molecular dynamics [32, 1, 50] can generate
time-resolved trajectories of protein folding and protein-ligand binding
predictions as well as structural ensemble simulations [54]. In these
simulations, mathematical models and numerical algorithms are imperative for
achieving computational accuracy and efficiency. A large number of advanced
algorithms have been developed to reduce the computational cost and improve
the accuracy for biomolecular simulations [52, 55, 10, 25]. A major difficulty
of molecular dynamics is the long timescales associated with real molecular
processes taking place in nature. Therefore, ignoring the requirement of
having time-resolved trajectories of the molecular processes will immediately
remove the difficulty. Indeed, it is sufficient for most studies to have a
predicted representative ensemble of structures for a given process. This
representative prediction can be generated by Monte Carlo sampling [21].

Monte Carlo method is one of the most of popular approaches for biomolecular
systems. Under a physiological condition, biomolecules are immersed in and
interact with surrounding water molecules and other possible co-factors. As
such, Monte Carlo simulations of a biomolecule have to deal with a large
number of solvent water molecules, which makes the simulations very expensive
and sometimes, intractable. Additionally, in Monte Carlo simulations, the
biomolecular conformation is subject to random perturbations [41]. These
perturbations will inevitably result in the overlaps between the biomolecule
and explicit solvent molecules, which leads to an unfavorable and non-
representative structure. Implicit solvent models, such as Poisson-Boltzmann
(PB) [15, 19], polarizable continuum [57, 14] and Generalized Born (GB)
methods [18, 49, 56, 42] are developed to overcome this challenge by taking a
mean field approximation of water molecules and resulting in a dielectric
continuum. The GB method is faster than PB methods but it only provides an
approximation for electrostatic energies. PB methods, derived from fundamental
physical theories [3, 48], offer more accurate electrostatic analysis. PB
model has been applied to the calculations of protein-protein and protein-
ligand binding energies[45], the pH value predictions of protonation and/or
deprotonation states of titration sites[30], and drug design [59]. To seek for
an accurate, efficient, and robust numerical solver, a large number of
numerical methods have been developed for the PB model, including finite
difference method (FDM) [29], finite element method (FEM) [2], and boundary
element method (BEM) [25, 39]. Among this variety of numerical explorations,
the FDM has the most enfranchisement such as Amber PBSA [61], Delphi [34],
APBS [2, 30], MIBPB [66, 23, 10, 46, 24], and CHARMM PBEQ [29]. Among them,
MIBPB is the solely available second-order accurate method and has been used
to calibrate the GB method in Amber [20], where PB methods are generally very
expensive. In addition, the molecular surface involved in all the
aforementioned method with corresponding software developed, such as ESES[36],
Nanoshaper[16], and MSMS[53].

Over the past a few years, machine learning, including deep learning, has had
tremendous success in science and engineering. Especially, convolutional
neural networks have proved their ability to automatically extract features
and recognize patterns from relatively simple but large datasets. Deep
learning has a growing dominance in important applications such as handwriting
recognition, speech recognition, and drug discovery [26, 40, 44]. Aided by the
availability of quality databases, new algorithms, graphics processing unit
(GPU), and high-performance computers, various machine learning approaches
have established in many classical computational problems such as solvation
free energies, protein-ligand binding affinities, mutation impacts, toxicity,
partition coefficients, protein B-factors, etc [33, 31, 28, 58, 47, 7, 6, 64,
65, 13]. Additionally, deep learning neural networks are also applied in
computational protein design [60], stability changes of protein induced by
mutations [8, 5], and calculations of protein energy[62, 11].

Recently, we developed a Poisson-Boltzmann based machine learning (PBML)
model, which can compute the solvation free energy of macromolecules in the
solvent with the GB speed and the PB accuracy [12]. We assume that all of the
macromolecular electrostatic solvation free energies follow a probability
distribution, which can be sampled by the PB model. Our idea is based on a
representability hypothesis and a learning hypothesis. The representability
hypothesis states that the solvation free energy of a molecule can be
described by the features of atom interactions and their geometric relations
in the solvent. Thus, we can construct feature vectors to characterize the
molecular electrostatic distribution. In our learning hypothesis, we assume
that a machine learning model can be trained based on training labels and
corresponding features for a sufficiently large training set of molecules.
Additionally, advanced machine learning algorithms can give accurate
predictions of the electrostatic potential for a new molecule which has the
same probability distribution with the training set. In our approach, training
labels are computed from MIBPB and features are generated using multiscale
weighted colored subgraphs [47].

In the present work, we apply our newly developed PBML model to compute
molecular solvation free energies in the implicit-solvent Monte Carlo
simulations, which typically require millions of samplings. The new machine
learning-based implicit-solvent Monte Carlo model can guarantee the accuracy
of the implicit-solvent Monte Carlo model while dramatically speeding up
existing implicit-solvent Monte Carlo algorithms.

This manuscript is organized as follows. Section II gives a brief introduction
of molecular force fields, Monte Carlo methods, and implicit solvent models.
The PBML model is introduced in this section as well, which includes the
Poisson-Boltzmann equation, Generalized Born model, and multiscale weighted
colored subgraphs. Section III presents the results of structural predictions
of benzene and the human hyperplastic discs protein (PDB: 1i2t) [17] in water.
We demonstrate that the PBML model is more accurate and faster than commonly
used PB solvers and thus, can significantly reduce the computational time of
implicit-solvent Monte Carlo simulations. A summary is given in Section IV.

##  II Methods and algorithms

In this section, we briefly review biomolecular force fields, the Monte Carlo
methods, and implicit solvent models, followed by the Poisson-Boltzmann based
machine learning model.

###  II.A Biomolecular force fields

The quality of molecular simulations depends crucially on molecular force
fields to offer a physical representation of molecular interactions and energy
distributions. Molecular force fields typically describe molecular
interactions in terms of classical molecular mechanics of atoms. The potential
energies of atomic interactions are approximated by a set of mathematical
functions, modeling the bonded and non-bonded components. These functions
consist of a set of free coefficients, which are obtained by approximating
either the results of elaborate quantum mechanical calculations, or
experimental data. One of the advantages of biomolecular force field approach
is its computational efficiency. The potential energy can be efficiently
computed at the molecular level comparing to other methods, such as quantum
mechanical approaches, which deal with electrons [9, 35]. Additionally, the
forces in molecular dynamics can be evaluated analytically from molecular
force fields.

A variety of molecular force fields have been developed for various purpose.
In this work, we adopt the popular and simple Amber ff99SB force field [35].
The Amber force field for governing the potential energy consists of the
following terms,

$E=\sum_{\rm bonds}k_{b}(r-r_{0})^{2}+\sum_{\rm
angles}k_{\theta}(\theta-\theta_{0})^{2}+\sum_{\rm
dihedrals}V_{n}[1+\cos(n\phi-\gamma)]+\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\Bigg{[}\frac{A_{ij}}{R_{ij}^{12}}-\frac{B_{ij}}{R_{ij}^{6}}+\frac{q_{i}q_{j}}{\epsilon_{1}R_{ij}}\Bigg{]}$
(1)

where $k_{b}$, $k_{\theta}$, and $V_{n}$ are force constants. Here, $r$,
$\theta$, and $\phi$ are bond length, angle and dihedral angle with $r_{0}$,
$\theta_{0}$, and $\gamma$ being optimal bond length, optimal angle and proper
dihedral angle, respectively. The first three terms in the energy expression
describe the bonded energy of the molecular system. The last term represents
the Lennard-Jones interactions and electrostatic interactions, where $N$ is
the number of atoms in the molecular system, $R_{ij}$ is the distance between
$i$th and $j$th atoms, $A_{ij}$ and $B_{ij}$ are Lennard-Jones parameters,
$q_{i}$ is the atom charge, and $\epsilon_{1}$ is the dielectric constant.

###  II.B Monte Carlo methods

In this session, we provide a brief introduction of the molecular dynamics and
the Monte Carlo method. We start from statistical mechanics and show that the
calculation of the physical property of a solute-solvent system using
molecular dynamics is computationally expensive or even intractable [21].
Then, we introduce Metropolis’s Monte Carlo method for biomolecular
simulations [41].

The classical expression for the partition function $Q$ of a solute-solvent
system is

$Q=c\int
d\textbf{r}d\textbf{p}\exp[-\mathcal{H}(\textbf{r},\textbf{p})/k_{B}T],$ (2)

where $\textbf{r}=\\{\textbf{X},\textbf{Y}\\}$ stands for the atomic
coordinates of a solute X and solvent Y, p stands for the corresponding
momenta, $c$ is a physical constant as specified below, $k_{B}$ is the
Boltzmann constant and $T$ is the temperature of the system. The function
$\mathcal{H}(\textbf{r},\textbf{p})$ is the Hamiltonian of the system. It
describes the total energy of an individual system as summation of the kinetic
energy $\mathcal{K}$ and the potential energy $E$:
$\mathcal{H}=\mathcal{K}+E$, where $\mathcal{K}$ is a quadratic function of
the momenta. For a system of $N$ identical atoms, one has $c=1/(h^{3N}N!)$
using the Planck constant $h$. Under the assumption that all of the other
physical observables $A$ of interest depend only on the positions, i.e.,
$A=A(\textbf{r})$, the integration over the momenta can be carried out
analytically in a classical mechanical treatment. As a result, the expected
value of a physical observable of interest is given by

$\langle A\rangle=\frac{\int d\textbf{r}A(\textbf{r})\exp[-\beta
E(\textbf{r})]}{\int d\textbf{r}\exp[-\beta E(\textbf{r})]},$ (3)

where $\beta=1/k_{B}T$. Evaluating $\langle A\rangle$ requires numerical
techniques, such as quadrature rules for the integration. Since each particle
moves in a three dimensional (3D) space, the total number of degrees of
freedom is $3N$ for a system of $N$ atoms. If each dimension is integrated
with a mesh size of $m$ points, the total number of points for the integration
is $m^{3N}$, which is computationally prohibitive.

The complexity in evaluating Eq. (3) can be significantly reduced by using the
Monte Carlo sampling. Indeed, Metropolis et al. [41] suggested an efficient
Monte Carlo scheme to approximate the ratio in Eq. (3). Let us denote the
probability density function in finding a microstate in the canonical ensemble
in a configuration r by

$P(\textbf{r})=\frac{\exp[-\beta E(\textbf{r})]}{\int d\textbf{r}\exp[-\beta
E(\textbf{r})]}.$ (4)

According to this probability function, we can perturb randomly selected
points in the configuration. Hence, the number of points $n_{i}$ generated per
unit volume in the neighborhood of r is equal to $N_{mc}\times P(\textbf{r})$
for the average of $A(\textbf{r})$, which is

$\langle
A\rangle\approx\frac{1}{N_{mc}}\sum_{i=1}^{N_{mc}}n_{i}A(\textbf{r}_{i}),$ (5)

where $N_{mc}$ is the total number running in Monte Carlo simulations.
Equation (5) shows that all states of ensemble contribute to the average
equally. Therefore, Metropolis Monte Carlo method starts at a given
configuration $\textbf{r}_{0}=\\{\textbf{X}_{0},\textbf{Y}_{0}\\}$ and next
perturbs the configuration by a defined transformation with a new
configuration $\textbf{r}_{1}=\\{\textbf{X}_{1},\textbf{Y}_{1}\\}$. The
probability to accept the new configuration is

$p_{acc}=\min(1,\exp(-\beta(E(\textbf{r}_{0})-E(\textbf{r}_{1})))).$ (6)

If the new configuration is rejected, the previous configuration is retained
and the method repeats another random perturbation. This process iterates
until the iteration number equals to a fixed number. It is shown that the
structure in the system will approach the Boltzmann distribution, if the
perturbations satisfy the condition

$\pi(\textbf{r}_{i})p_{ij}=\pi(\textbf{r}_{j})p_{ji},$ (7)

where $\pi(\textbf{r}_{i})$ is the probability of the system in configuration
$\textbf{r}_{i}$ and $p_{ij}$ is the probability to perturb the configuration
from state $\textbf{r}_{i}$ to state $\textbf{r}_{j}$ [41].

###  II.C Implicit solvent models

Implicit solvent models are class of multiscale techniques for reducing the
dimensionality of a solvent-solute system. They retain the crucial
electrostatic interactions between a biomolecule and its solvent environment
without modeling solvent molecules explicitly. A variety of two-scale implicit
solvent models has been developed, such as the Poisson-Boltzmann (PB) model
[19] and the generalized Born (GB) model [18, 49, 56, 42]. One desirable
application of implicit solvent models is the Monte Carlo simulations of
biomolecule in solvent, which is relatively easy to implement. The basic
derivation for molecular implicit solvent models relies on statistical
mechanics. For more detail, the reader is referred to the literature [51].
Essentially, the molecular solvation free energy can be given by

$\Delta G_{\rm solv}=\Delta G_{\rm elec}+\Delta G_{\rm nonpol},$ (8)

where $\Delta G_{\rm elec}$ represents the electrostatic contribution of the
solvent-solute interaction, and $\Delta G_{\rm nonpol}$ denotes the nonpolar
energy in the reversible work needed to insert a fixed configuration molecule
into the solvent with all solute charges set to zero. Here $\Delta G_{\rm
nonpol}$ is proportional to the solvent accessible surface area. The molecular
solvation free energy is used in our implicit-solvent Monte Carlo method to
represent solvent-solute interactions.

###  II.D  Poisson-Boltzmann based machine learning (PBML) model

In this section, we briefly discuss the Poisson-Boltzmann based machine
learning (PBML) model [12], which is applied to compute $\Delta G_{\rm elec}$
in Eq. (8). Our PBML model involves three major components, i.e., training
labels, molecular features, and learning algorithms. Our training labels for a
large training set of molecules are generated from solving the Poisson-
Boltzmann (PB) equation. Our molecular features for both the training set and
the test set constitute two parts, a GB part and a correction part. The latter
is computed from multiscale weighted colored subgraphs [12].

#### The Poisson-Boltzmann (PB) model

The PB model considers the solute biomolecule with $N_{c}$ fixed charges as
the interior domain $\Omega_{1}$, and the solvent, including free ions, as the
exterior domain $\Omega_{2}$. The interface $\Gamma$ separates these two
domains. The PB model is given as

$-\nabla\cdot\epsilon({\bf r})\nabla\phi({\bf r})+\bar{\kappa}^{2}({\bf
r})\phi({\bf r})=\sum_{k=1}^{N_{c}}q_{k}\delta({\bf r}-{\bf r}_{k}),$ (9)

For ${\bf r}\in\mathbb{R}^{3}$, $\phi(\bf{r})$ is the electrostatic potential,
$\epsilon(\bf{r})$ dielectric constant given by

$\epsilon(\bf{r})=\left\\{\begin{array}[]{ll}\epsilon_{1},\quad{\bf{r}}\in\Omega_{1},\\\
\epsilon_{2},\quad{\bf{r}}\in\Omega_{2},\end{array}\right.$ (10)

In the PB model, $\bar{\kappa}$ is the screening parameter with the relation
$\bar{\kappa}^{2}=\epsilon_{2}\kappa^{2}$ where $\kappa$ is the inverse Debye
length measuring the ionic effective length. To ensure the continuity of
electrostatic potential and flux density across the interface $\Gamma$, the PB
equation is associated with following interface conditions

$\phi_{1}({\bf r})=\phi_{2}({\bf
r}),\quad\epsilon_{1}\frac{\partial\phi_{1}({\bf r})}{\partial{\bf
n}}=\epsilon_{2}\frac{\partial\phi_{2}({\bf r})}{\partial{\bf n}},\quad{\bf
r}\in\Gamma$ (11)

where $\phi_{1}$ and $\phi_{2}$ are electrostatic potential from the solute
domain $\Omega_{1}$ and the solvent domain $\Omega_{2}$, and $\bf{n}$ is the
outward unit normal vector on $\Gamma$.

The solvation free energy can be obtained from the PB model by

$\Delta G_{\rm elec}^{\rm
PB}=\frac{1}{2}\sum_{k=1}^{N_{c}}q_{k}(\phi(\bf{r}_{k})-\phi_{0}(\bf{r}_{k}))$
(12)

where $\phi_{0}(\bf{r}_{k})$ is the free space solution to the PB equation
assuming no solvent-solute interface. To solve the PB equation, we apply the
accurate and robust 2nd order MIBPB solver [10, 46] developed in our group,
which applies rigorous treatment on geometric complexity, interface condition,
and charge singularity. The $\Delta G_{\rm elec}^{\rm PB}$ results generated
by MIBPB solver for a set of macromolecules are used as the training labels in
the representability hypothesis.

#### The Generalized Born (GB) model

Having described the labels for our machine learning training, we discuss the
molecular feature construction for both machine learning training and test,
which involves the GB model. As a fast approximation to the PB model, the GB
model compute the electrostatic solvation free energy by

$\Delta G^{\rm GB}_{\rm elec}\approx\sum_{i,j}\Delta G^{\rm
GB}_{ij}=-\frac{1}{2}\Big{(}\frac{1}{\epsilon_{1}}-\frac{1}{\epsilon_{2}}\Big{)}\frac{1}{1+\alpha\beta}\sum_{i,j}q_{i}q_{j}\Big{(}\frac{1}{f_{ij}(r_{ij},R_{i},R_{j})}+\frac{\alpha\beta}{B}\Big{)},$
(13)

where $R_{i}$ is the effective Born radius for $i$-th atom, $r_{ij}$ is the
distance between atoms $i$ and $j$, $\beta=\epsilon_{1}/\epsilon_{2}$,
$\alpha=0.571412$, and $B$ is the electrostatic size of the molecule. The
function $f_{ij}$ is given as

$f_{ij}=\sqrt{r^{2}_{ij}+R_{i}R_{j}{\rm
exp}\Big{(}-\frac{r^{2}_{ij}}{4R_{i}R_{j}}\Big{)}}.$ (14)

The effective Born radii $R_{i}$ is calculated by the following boundary
integral

$R^{-1}_{i}={\Big{(}-\frac{1}{4\pi}\oint_{\Gamma}\frac{{\bf r}-{\bf
r}_{i}}{|{\bf r}-{\bf r}_{i}|^{6}}\cdot\text{d}\bf{S}\Big{)}}^{1/3}.$ (15)

In Eq. (15), the MSMS package [43] is used to generate the triangulation
discretization of the molecular surface for the numerical surface integral on
$\Gamma$.

#### Multiscale weighted colored subgraphs

The weighted colored subgraph (WCS) use the notion $G(V,E)$ with vertices $V$
and edges $E$ to describe the atomic interactions in a protein of $N$ atoms.
The vertices is defined as

$\displaystyle V=\\{({\bf r}_{i},\alpha_{i})|{\bf
r}_{i}\in\mathbb{R}^{3},\alpha_{i}\in\mathcal{C},i=1,2,\dots,N\\},$ (16)

where $\mathcal{C}=\\{{\rm C,N,O,S,H}\\}$ contains all the commonly occurring
element types in a protein. Each vertex is an atom labeled by both its
position ${\bf r}_{i}$ element type $\alpha_{i}$, for $i=1,\cdots N$.

The edge $E$ relates the pairwise interactions, which are defined as a colored
set $\mathcal{P}=\\{\alpha\beta\\}$ with $\alpha,\beta\in\mathcal{C}$. For
$\mathcal{C}$ defined above, $\mathcal{P}=\\{{\rm
CC,CN,CO,CS,CH,NN,NO,NS,NH,OO,OS,OH,SS,SH,HH}\\}$ and we defined the partition
of $\mathcal{P}$ as $\mathcal{P}_{k}$, $k=1,2,...,15$ such that
$\mathcal{P}_{1}=\\{{\rm CC}\\}$, $\mathcal{P}_{2}=\\{{\rm CN}\\}$ and so on.
The set of involved vertices $V_{{\cal P}_{k}}$ is a subset of $V$ containing
all atoms involved in forming the pair in $\mathcal{P}_{k}$. For instance,
$\mathcal{P}_{2}=\\{{\rm CN}\\}$ contains all carbon-nitrogen atom pairs and
$V_{{\cal P}_{2}}$ contains all carbon and nitrogen atom vertices in the
protein. Based on these configuration, all the edges for pairwise atomic
interactions in the WCS description are defined by

$\displaystyle
E_{\mathcal{P}_{k}}^{\sigma,\tau,\zeta}=\\{\Phi^{\sigma}_{\tau,\zeta}(\|{\bf
r}_{i}-{\bf
r}_{j}\|)~{}|~{}\alpha_{i}\beta_{j}\in\mathcal{P}_{k};i=1,2,\dots,N_{\alpha},j=1,2,\dots,N_{\beta}\\},$
(17)

where $\|{\bf r}_{i}-{\bf r}_{j}\|$ defines the Euclidean distance between
$i^{th}$ and $j^{th}$ atoms, $N_{\alpha}$ and $N_{\beta}$ are numbers of type
$\alpha$ and $\beta$ atoms, $\sigma$ indicates the type of radial basic
functions (e.g., $\sigma=\text{L}$ for Lorentz kernel, $\sigma=\text{E}$ for
exponential kernel), $\tau$ is a scale distance factor between two atoms and
$\zeta$ is a parameter of power in the kernel (i.e., $\zeta=\kappa$ for
$\sigma=\text{E}$, $\zeta=\nu$ for $\sigma=\text{L}$). In this model, we use
generalized exponential functions

$\displaystyle\Phi^{\rm
E}_{\tau,\kappa}=e^{-(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|/\tau(r_{i}+r_{j}))^{\kappa}},\quad\kappa>0,$
(18)

and generalized Lorentz functions

$\displaystyle\Phi^{\rm
L}_{\tau,\nu}(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|)=\dfrac{1}{1+(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|/\tau(r_{i}+r_{j}))^{\nu}},\quad\nu>0,$
(19)

where $r_{i}$ and $r_{j}$ are, respectively, the van der Waals radius of the
$i^{th}$ and $j^{th}$ atoms. Finally, the features for describe the
electrostatics interactions and geometric properties are expressed as

$\mu^{k,\sigma,\tau,\nu,w}=\sum_{i=1}^{N_{\alpha}}\sum_{j=1}^{N_{\beta}}w_{ij}\Phi^{\sigma}_{\tau,\nu}(\|\mathbf{r}_{i}-\mathbf{r}_{j}\|),\quad\alpha_{i}\beta_{j}\in\mathcal{P}_{k},$
(20)

where $w_{ij}$ is a weight function assigned to each atomic pair with
$w_{ij}=1$ for atomic rigidity or $w_{ij}=q_{j}$ for atomic charge. Since we
have 15 options of the colored subsets $\mathcal{P}_{k}$, we can obtain
corresponding 15 subgraph centralities $\mu^{k,\sigma,\tau,\nu,w}$, for
$k=1,2,\dots,15$. By varying kernel parameters $(\sigma,\tau,\nu,w)$, one can
achieve multiscale centralities for multiscale weighted colored subgraph
(MWCS) [4], which can be the features.

With labels and features described above, we can construct the machine
learning model to predict the solvation free energy of new macromolecules.
Specifically, using MIBPB results as labels, and GB and MWCS results as
features, we train gradient boosting decision trees (GBDTs) for the solvation
free energy prediction.

##  III Results

In this section, we demonstrate the performance of the proposed MLIMC method
numerically. First, we describe the Poisson-Boltzmann based machine learning
(PBML) model for computing protein electrostatic solvation energies, followed
by the illustration of the accuracy and efficiency of the model. The use of
the PBML model for electrostatic interactions in the MC simulations is
introduced. Our main idea is to replace time-consuming electrostatic
calculations by using our PBML model. The efficiency of our new MLIMC model is
also examined. Finally, we validate the proposed MLIMC method by two cases.
Case 1 is a small molecule, benzene, with initial atom position randomly
protruded. Our MLIMC method is used to reconstruct the benzene molecule in
solvent. Case 2 is a relatively larger molecule, protein (PDB: 1i2t) with 61
amino acid residues. In this case, we stretch the last two residues of 1i2t
using steered molecular dynamics and then we try to restore the equilibrium
configuration by using the proposed MLIMC method. Both simulations are carried
out at a temperature of $27^{\circ}\text{C}$, the dielectric constants are
$\epsilon_{1}=1$ in the molecule and $\epsilon_{2}=80$ in the solvent, the
MSMS [43] mesh density is set as 2, and the Debye-Huckel constant is set as
$\kappa=0.1257$Å-1. There are three kernels used to generated features for
machine learning, which are $(\text{E},0.3,2,1)$, $(\text{E},4.7,2,q_{j})$,
and $(\text{L},4.2,5,1)$.

To measure the performance, we use the root-mean-square deviation (RMSD) of
atomic positions in length units (Å), defined as

$\text{RMSD}(\textbf{v},\textbf{w})=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\Big{(}(v_{ix}-w_{ix})^{2}+(v_{iy}-w_{iy})^{2}+(v_{iz}-w_{iz})^{2}\Big{)}},$
(21)

where $\textbf{v},\textbf{w}\in\mathbb{R}^{N\times 3}$ are vectors of
positions of the $N$ atoms at two different MC samplings. Moreover, we also
present relative errors of the total energy measured by comparing the energy
for a MC sampling $E_{MC}$, and the energy for the equilibrium state $E_{SS}$
as

$\text{e}_{\text{ttl}}=\frac{|E_{SS}-E_{MC}|}{|E_{SS}|}\times 100\%.$ (22)

We compute the RMSD and errors between Monte Carlo sampling results and the
original molecular structure for every $100$ Monte Carlo steps for both cases.
The core code was written in C/C++ and a cython wrapper calling the core code
for performing adds-on functions and applications. Our simulations are
produced on a desktop with an i5 7500 CPU and 16GB memory.

###  III.A PBML model

The MLPB model used in Monte Carlo simulation is a pre-trained model. The
training set includes 3706 protein structures from the PDBbind v2015 refined
set[37]. This refined set was selected from a general set of 14,620 protein-
ligand complexes. A data pre-processing (i.e. adding force field parameters)
is required before a PB solver can be used for electrostatics calculations.
Though the PDBbind refined set consists of protein-ligand complexes, only
protein structures are applied for calculations. These protein structures are
adjusted by the protein preparation wizard utility of the Schrodinger 2015-2
Suite [38] with default parameters unless filling the missing side chains is
required.

The training set covers a wide range of proteins in different sizes with atom
numbers from 997 to 27,713. The current training set can be expanded to an
even larger group of proteins. However, from our test, we conclude that
expanding training set will not significantly improve the trained model, thus
the size of the current training set is sufficiently large.

The purpose of PBML is to implement a machine learning predictor of PB
electrostatic solvation free energies for various proteins efficiently and
accurately without explicitly solving the PB equation. Gradient boosting
decision tree method is selected for this supervised learning task because of
its efficiency. To accuracy of the PBML model is maintained by the accurate
electrostatic free energy of solvation as the label calculated by the MIBPB
solver. Once a trained PBML model is obtained, the MIBPB solver will not be
called anymore. Using the learned PBML model only requires calculating
features on the prediction of electrostatic solvation free energies for new
compounds, which is rapid.

###  III.B Efficiency of the PBML model

Figure 1:  Comparison of the mean CPU time (second) per protein and the mean
absolute percentage errors of Amber, DelPhi and machine learning predictions
of the electrostatic solvation free energies using the test set of 195
proteins. Left: Results of Amber and DelPhi were obtained at ten different
mesh sizes from 0.2 to 1.1 Å; Results of PBML were obtained at four MSMS
densities (number of vertices per Å2) at 15, 2, 1, and 0.5. Right: A zoom-in
plot of the left plot for small CPU time.

Figure 1 shows the results for computing solvation energy on 195 proteins from
PDBbind v2015 core set [37] using PBML, Amber, and Dephi. The results are
shown in terms of the average CPU time per protein versus the mean absolute
percentage errors. From Fig. 1(a), we can see PBML is more accurate and much
faster than standard PB solvers such as DelPhi and Amber PB. Figure 1(b) gives
more details by zooming into the region where CPU time is small to distinct
the CPU time used by the PBML using different MSMS density.

We here add a few notes about how we improve the PBML model in addition to
machine learning. We notice that in the energy and feature calculations, every
term has a degree of freedom associated with the number of atoms, except the
computation of the effective Born radii $R_{i}$ in Eq. (15), which depends on
the number of surface triangles $M$. Since $M\gg N$, faster evaluating of Eq.
(15) can significantly accelerate the entire Monte Carlo process. In our
present implementation, instead of taking the integral in Eq (15) on each
triangle, we take the integral on a neighborhood of each vertex. This
treatment nearly doubled the efficiency of the GB method since number of
vertices is about half of number of triangles on the surface. In addition,
applying a cut-off can also further improve the GB method.

###  III.C MLIMC model

The assembling of MLIMC includes the implementation of empirical potential
energy functions (except electrostatics) and the prediction of electrostatics
for each step on Monte Carlo simulations. The conformation of the target
protein is perturbed randomly on each step. The new conformation is directly
accepted if it shows a lower energy or is accepted with a probability
determined by the Boltzmann distribution if it shows a higher energy. As the
MLPB model is pre-trained before simulations, the Monte Carlo simulation does
not include the time for solving the PB equation, resulting in much reduced
time for MLIMC simulations.

###  III.D Efficiency of the MLIMC model

We show that the high efficiency of the MLPB model will significantly improve
the efficiency of the MLIMC model.

Table 1: Average CPU time for one step MLIMC simulation using Amber, DelPhi and PBML for electrostatic solvation free energy on the 195 protein dataset. Results of Amber and DelPhi were obtained at $0.2\text{\AA}$ and $0.5\text{\AA}$ mesh sizes, and that from PBML uses mesh density $2$. The average CPU time includes all computations needed for Monte Carlo evaluations. The PB error is obtained relative to the electrostatic solvation energy computed from MIBPB solver with grid size $h=0.2\text{\AA}$.  |  | $h=0.2\text{\AA}$ | $h=0.5\text{\AA}$
---|---|---|---
PB solvers | PBML | Amber | DelPhi | Amber | DelPhi
Avg. CPU time (sec) | 25 | 6136 | 1621 | 1177 | 214
PB error (%) | 0.484 | 0.618 | 0.819 | 1.271 | 1.552

Table 1 shows the mean CPU time of one Monte Carlo step and the mean absolute
percentage errors of Amber, DelPhi and PBML predictions of the electrostatic
solvation free energies of the 195 proteins. The mean CPU time for each
protein includes the computations for the total energies, in which computing
electrostatic is the dominant component. Clearly, the machine learning method
has the highest accuracy but lowest CPU time. For the same accuracy level
($<1\%$), the estimation of the mean CPU time for a one-million-step Monte
Carlo simulation is $6.136\times 10^{8}$s, $1.621\times 10^{8}$s and
$2.5\times 10^{6}$s for using Amber, DelPhi and PBML, respectively. Even with
compromised accuracy for DelPhi and Amber at gird size of $0.5\text{\AA}$, the
MLIMC with PBML will be $47$ times faster than that with Amber and $8$ times
faster than that with DelPhi. Next we show some MC simulation results using
MLIMC on the benzene molecule and the human hyperplastic discs protein
(PDB:1i2t).

###  III.E Test case one: Benzene molecule

(a)(b)

Figure 2: MLIMC simulation of benzene in solvent. (a) the red curve is the
total energy calculated by our implicit-solvent Monte Carlo model and the blue
curve is the root mean square deviation of the atomic positions on each Monte
Carlo step to the non-protruded one. (b) the red curve is the error of total
energies e${}_{\text{ttl}}$ defined by Eq. 22 and the blue curve is the same
RMSD as the left figure.

(a)(b)

Figure 3: Illustration of MLIMC simulations of a benzene molecular in solvent.
(a) the blue structure is the randomly perturbed atom positions and the green
one is the benzene structure in steady state. (b) the red one is the benzene
structure after MLIMC sampling compared with the equilibrium structure in
green. Pictures are produced with VMD [27].

Our first case is a Benzene molecule with some atomic position randomly
perturbed. In detail, we fixed three atoms at equilibrium positions in order
to have the prediction and the comparison structure in the same plane, and
perturb the coordinates of the remained nine atoms in $(\rho,\theta,\phi)$
directions by uniformly distributed random numbers in
$([0,10],[0,2\pi],[0,\pi])$. The initial RMSD is $6.42\text{\AA}$ as compared
with the equilibrium position. We will try to perform a MC simulation on this
perturbed molecule to see if the original steady status can be obtained.
Figure 2(a) shows the total energy and RMSD vs MC steps, from which we can see
that the total energy of benzene in solvent starts at $349123.61$ kcal/mol and
converges to the range of $(5,15)$ kcal/mol after the first $20,000$ MC steps.
It stays in a convergent range for the rest MC steps. The RMSD initially is
$6.42\text{\AA}$ and ends around $0.15\text{\AA}$. It decreases rapidly as the
total energy for the first $20,000$ steps. After $20,000$ steps, the total
energy converges with only slightly oscillation, and the RMSD keeps the
decreasing trend until it reaches around $0.15\text{\AA}$ when MC steps are
greater than $70,000$.

Figure 2(b) shows errors and RMSD versus MC steps. Here we set $E_{SS}$ in Eq.
(22) to be $10.60$ kcal/mol as the steady state energy for reference. The plot
shows that the errors of total energy are very small for our MC simulation
after $10,000$ iterations. When the simulation structure is close to that of
its equilibrium state, the RMSD is smaller than $1\text{\AA}$ and the errors
stay in between $1\%$ to $100\%$. Note since the total energy is a small
number, a tiny perturbation causes a large error changing.

Qualitatively, Fig. 3(a) shows that the benzene molecule with its initial
perturbed structure in blue and the equilibrium structure is in green. After
the MC simulation, we receive the predicted structure in red as compared with
the steady state structure in green as shown in Fig. 3(b). The total CPU time
for $100,000$ Monte Carlo steps is $643$ seconds.

###  III.F Test case 2: Protein (PDB: 1i2t)

(a)(b)

Figure 4: MLIMC simulation of the protein (PDB: 1i2t) in solvent. (a) the red
curve is the total energy calculated by implicit-solvent Monte Carlo model,
the blue curve is the root mean square deviation of the atomic positions on
each Monte Carlo step to the non-protruded one; (b) the red curve is the error
of total energies e${}_{\text{ttl}}$ defined by Eq. 22, the blue curve is the
same RMSD as the left figure.

(a)(b)

Figure 5: Illustration of MLIMC simulations of the protein (PDB: 1i2t) [17] in
solvent. (a) the blue structure is the perturbed protein structure generated
with steered molecular dynamics and the green one is the original structure at
the equilibrium state in solvent. (b) the red one is the predicted structure
after MLIMC sampling compared with the original structure in green. Pictures
are produced with VMD [27].

The second MC test is on the human hyperplastic discs protein (PDB: 1i2t) with
61 residues. We first stretch the last two residues of the original protein by
a steered molecular dynamics. As a result, the stretched molecule has an
initial RMSD of $8.14\text{\AA}$. We apply our MLIMC for $100,000$ steps,
which takes $16,684$ seconds in CPU time. Figure 4(a) shows that the total
energy of $7260.90$ kcal/mol initially decays rapidly within the first $5,000$
Monte Carlo steps, then oscillates around $-2070.00$ kcal/mol. In the same
plot, we can see the RMSD drops quickly in the first $10,000$ MC steps, after
than decays slowly with fluctuation for the next $40,000$ MC steps, then
decays steadily after $45,000$ MC steps, and finally oscillates slightly
around $2.2\text{\AA}$ after $60,000$ MC steps. For the energy errors shown in
Fig. 4 (b), relative to the total energy in the equilibrium of $-2068.13$
kcal/mol, the errors rapidly decays in the first $20,000$ MC steps and then
oscillate within 10% after that.

Similar to the Benzene case, Fig. 5(a) qualitatively shows the perturbed
structure in blue against the steady state structure in green for the first
two residues and Fig. 5 (b) shows that the MLIMC structural prediction in red
color after the MC simulation, which is very close to the steady state
structure in green.

##  IV Conclusion

Monte Carlo simulations are widely used in science and engineering for
molecular structure optimization and prediction. In many situations,
particularly biomolecular systems, the solute molecule is immersed in a water
solvent and the full-scale explicit solvent Monte Carlo simulations are very
expensive. Alternatively, implicit solvent Monte Carlo methods using either
Poisson-Boltzmann (PB) model or generalized Born (GB) model for computing
electrostatics can greatly reduce the degree of freedom. However, the accuracy
reduction in GB model or the efficiency concerns in PB model hinders the wide
application of implicit solvent Monte Carlo simulation. In this work, we
introduce a machine learning-based implicit-solvent Monte Carlo (MLIMC) method
for molecular structure optimization and prediction. A vital component of our
MLIMC is the newly developed Poisson-Boltzmann based machine learning (PBML)
model, which maintains the PB accuracy at the GB cost. We validate the
proposed MLIMC method by simulating two molecular systems, randomly perturbed
benzene structure and protein (PDB: 1i2t) structures modified by a steered
molecular dynamics. Numerical experiments demonstrate that proposed MLIMC is
efficient in predicting molecular structures at equilibrium. In a comparative
analysis, we show that the MLIMC model has a great advantage on CPU time and
accuracy over DelPhi and Amber PB based Monte Carlo methods. We believe this
innovated PBML method can also disruptively change the current status of PB
based molecular simulation involving molecular dynamics [22] and Monte Carlo.
MLIMC provides accurate electrostatic solvation energy at each configuration
of the target protein thus can be helpful in searching protein folding states
as intermediate or final using MC based simulation. The resulting machine
learning-based implicit molecular dynamics (MLIMD), together with the present
MLIMC model, will have a vast variety of applications in molecular science,
including drug design.

## Acknowledgment

This work was supported in part by NIH grant GM126189, NSF grants DMS-2052983,
DMS-1761320, and IIS-1900473, NASA grant 80NSSC21M0023, Michigan Economic
Development Corporation, MSU Foundation, Bristol-Myers Squibb 65109, and
Pfizer. The work of WG was supported in part by NSF grants DMS-1819193 and
DMS-2110922.

## Literature cited

##

  * [1] B. J. Alder and T. E. Wainwright.  Studies in Molecular Dynamics. I. General Method.  J Chem Phys, 31:459–466, 1959.
  * [2] N. A. Baker, D. Sept, M. J. Holst, and J. A. Mccammon.  The adaptive multilevel finite element solution of the {{Poisson-Boltzmann}} equation on massively parallel computers.  IBM Journal of Research and Development, 45(3-4):427–438, 2001\.
  * [3] D. Beglov and B. Roux.  Solvation of complex molecules in a polar liquid: an integral equation theory.  Journal of Chemical Physics, 104(21):8678–8689, 1996.
  * [4] D. Bramer and G.-W. Wei.  Multiscale weighted colored graphs for protein flexibility and rigidity analysis.  The Journal of chemical physics, 140:054103, 2018.
  * [5] Z. Cang and G.-W. Wei.  Analysis and prediction of protein folding energy changes upon mutation by element specific persistent homology.  Bioinformatics, 33(22):3549–3557, 2017.
  * [6] Z. X. Cang, L. Mu, and G. W. Wei.  Representability of algebraic topology for biomolecules in machine learning based scoring and virtual screening .  PLOS Computational Biology, 14(1):e1005929, https://doi.org/10.1371/journal.pcbi.100, 2018.
  * [7] Z. X. Cang and G. W. Wei.  Integration of element specific persistent homology and machine learning for protein-ligand binding affinity prediction .  International Journal for Numerical Methods in Biomedical Engineering, 34(2):DOI: 10.1002/cnm.2914, 2018.
  * [8] H. Cao, J. Wang, L. He, Y. Qi, and J. Z. Zhang.  Deepddg: predicting the stability change of protein point mutations using neural networks.  Journal of chemical information and modeling, 59(4):1508–1514, 2019\.
  * [9] A. Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev and R. J. Simmerling, C.; Wang, B.; Woods.  The Amber Biomolecular Simulation Programs.  J. Comput. Chem., 26:1668–1688, 2005.
  * [10] D. Chen, Z. Chen, C. Chen, W. H. Geng, and G. W. Wei.  MIBPB: A software package for electrostatic analysis.  J. Comput. Chem., 32:657–670, 2011.
  * [11] J. Chen, X. Xu, S. Liu, and D. H. Zhang.  A neural network potential energy surface for the f+ ch 4 reaction including multiple channels based on coupled cluster theory.  Physical Chemistry Chemical Physics, 20(14):9090–9100, 2018.
  * [12] J. Chen, Y. Xu, Z. Cang, W. Geng, and G.-W. Wei.  Poisson-Boltzmann based machine learning (PBML) model for electrostatic analysis.  Preprint, 2021.
  * [13] T. Ching, D. S. Himmelstein, B. K. Beaulieu-Jones, A. A. Kalinin, B. T. Do, G. P. Way, E. Ferrero, P.-M. Agapow, M. Zietz, M. M. Hoffman, and Others.  Opportunities and obstacles for deep learning in biology and medicine.  Journal of The Royal Society Interface, 15(141):20170387, 2018.
  * [14] M. Cossi, V. Barone, R. Cammi, and J. Tomasi.  Ab initio study of solvated molecules: A new implementation of the polarizable continuum model.  Chemical Physics Letters, 255:327–335, 1996.
  * [15] M. E. Davis and J. A. McCammon.  Electrostatics in biomolecular structure and dynamics.  Chemical Reviews, 94:509–521, 1990.
  * [16] S. Decherchi and W. Rocchia.  A general and robust ray-casting-based algorithm for triangulating surfaces at the nanoscale.  PloS one, 8(4):e59744, 2013.
  * [17] R. C. Deo, N. Sonenberg, and S. K. Burley.  X-ray structure of the human hyperplastic discs protein: an ortholog of the c-terminal domain of poly (a)-binding protein.  Proceedings of the National Academy of Sciences, 98(8):4414–4419, 2001.
  * [18] B. N. Dominy and C. L. Brooks III.  Development of a generalized Born model parameterization for proteins and nucleic acids.  Journal of Physical Chemistry B, 103(18):3765–3773, 1999.
  * [19] F. Fogolari, A. Brigo, and H. Molinari.  The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology.  Journal of Molecular Recognition, 15(6):377–392, 2002.
  * [20] N. Forouzesh, S. Izadi, and A. V. Onufriev.  Grid-based surface generalized born model for calculation of electrostatic binding free energies.  Journal of chemical information and modeling, 57(10):2505–2513, 2017.
  * [21] D. Frenkel.  Introduction to Monte Carlo Methods.  NIC Series, 23:29–59, 2004.
  * [22] W. Geng and G.-W. Wei.  Multiscale molecular dynamics using the matched interface and boundary method.  Journal of computational physics, 230(2):435–457, 2011.
  * [23] W. Geng, S. Yu, and G. W. Wei.  Treatment of charge singularities in implicit solvent models.  Journal of Chemical Physics, 127:114106, 2007.
  * [24] W. Geng and S. Zhao.  A two-component matched interface and boundary (mib) regularization for charge singularity in implicit solvation.  Journal of Computational Physics, 351:25–39, 2017.
  * [25] W. H. Geng and R. Krasny.  A treecode-accelerated boundary integral Poisson-Boltzmann solver for continuum electrostatics of solvated biomolecules.  J. Comput. Phys., 247:62–87, 2013.
  * [26] T. B. Hughes, G. P. Miller, and S. J. Swamidass.  Modeling epoxidation of drug-like molecules with a deep machine learning network.  2015\.
  * [27] W. Humphrey, A. Dalke, and K. Schulten.  VMD – Visual Molecular Dynamics.  Journal of Molecular Graphics, 14(1):33–38, 1996.
  * [28] J. Jiménez, M. Skalic, G. Mart$\backslash$’$\backslash$inez-Rosell, and G. De Fabritiis.  K DEEP: Protein–Ligand Absolute Binding Affinity Prediction via 3D-Convolutional Neural Networks.  Journal of chemical information and modeling, 58(2):287–296, 2018\.
  * [29] S. Jo, M. Vargyas, J. Vasko-Szedlar, B. Roux, and W. Im.  PBEQ-Solver for online visualization of electrostatic potential of biomolecules.  Nucleic Acids Research, 36:W270 –W275, 2008.
  * [30] E. Jurrus, D. Engel, K. Star, K. Monson, J. Brandi, L. Felberg, D. Brookes, L. Wilson, J. Chen, K. Liles, M. Chun, P. Li, D. Gohara, T. Dolinsky, R. Konecny, D. Koes, J. Nielsen, T. Head-Gordon, W. Geng, R. Krasny, G.-W. Wei, M. Holst, J. McCammon, and N. Baker.  Improvements to the APBS biomolecular solvation software suite.  Protein Science, 27(1), 2018.
  * [31] M. Karimi, D. Wu, Z. Wang, and Y. Shen.  DeepAffinity: Interpretable Deep Learning of Compound-Protein Affinity through Unified Recurrent and Convolutional Neural Networks.  arXiv preprint arXiv:1806.07537, 2018.
  * [32] J. Karplus, M.; Kuriyan.  Molecular Dynamics and Protein Function.  Proc. Nat. Acad. Sci. USA, 102:6679–6685, 2005.
  * [33] A. Korotcov, V. Tkachenko, D. P. Russo, and S. Ekins.  Comparison of deep learning with multiple machine learning methods and metrics using diverse drug discovery data sets.  Molecular pharmaceutics, 14(12):4462–4475, 2017.
  * [34] L. Li, C. Li, S. Sarkar, J. Zhang, S. Witham, Z. Zhang, L. Wang, N. Smith, M. Petukh, and E. Alexov.  DelPhi: a comprehensive suite for DelPhi software and associated resources.  BMC biophysics, 5(1):9, 2012.
  * [35] D. E. Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw.  Improved Side-Chain Torsion Potentials for the Amber ff99SB Protein Force Field.  Proteins, 78:1950–1958, 2010.
  * [36] B. Liu, B. Wang, R. Zhao, Y. Tong, and G.-W. Wei.  Eses: Software for e ulerian solvent excluded surface, 2017.
  * [37] Z. Liu, Y. Li, L. Han, J. Liu, Z. Zhao, W. Nie, Y. Liu, and R. Wang.  PDB-wide collection of binding data: current status of the PDBbind database.  Bioinformatics, 31(3):405–412, 2015.
  * [38] S. LLC.  Schrödinger release 2015-2, Schrödinger LLC, new york, 2015\.
  * [39] B. Lu, X. Cheng, J. Huang, and J. A. McCammon.  AFMPB: An Adaptive Fast Multipole Poisson-Boltzmann Solver for Calculating Electrostatics in Biomolecular Systems.  Comput. Phys. Commun., 184:2618–2619, 2013.
  * [40] A. Lusci, G. Pollastri, and P. Baldi.  Deep architectures and deep learning in chemoinformatics: the prediction of aqueous solubility for drug-like molecules.  Journal of chemical information and modeling, 53(7):1563–1575, 2013\.
  * [41] E. Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller.  Equation of State Calculations by Fast Computing Machines.  J. Chem. Phys., 21:1087–1092, 1953.
  * [42] J. Mongan, C. Simmerling, J. A. McCammon, D. A. Case, and A. Onufriev.  Generalized Born Model with a Simple, Robust Molecular Volume Correction.  Journal of Chemical Theory and Computation, 3(1):159–169, 2007\.
  * [43] MSMS.  https://mgl.scripps.edu/people/sanner/html/msms_home.html.
  * [44] D. D. Nguyen, Z. Cang, K. Wu, M. Wang, Y. Cao, and G.-W. Wei.  Mathematical deep learning for pose and binding affinity prediction and ranking in d3r grand challenges.  Journal of computer-aided molecular design, 33(1):71–82, 2019.
  * [45] D. D. Nguyen, B. Wang, and G.-W. Wei.  Accurate, robust, and reliable calculations of poisson–boltzmann binding energies.  Journal of computational chemistry, 38(13):941–948, 2017.
  * [46] D. D. Nguyen, B. Wang, and G. W. Wei.  Accurate, robust and reliable calculations of {Poisson-Boltzmann} binding energies.  Journal of Computational Chemistry, 38:941–948, 2017.
  * [47] D. D. Nguyen, T. Xiao, M. L. Wang, and G. W. Wei.  Rigidity strengthening: A mechanism for protein-ligand binding .  Journal of Chemical Information and Modeling, 57:1715–1721, 2017\.
  * [48] A. Onufriev, D. Bashford, and D. A. Case.  Modification of the generalized {Born} model suitable for macromolecules.  Journal of Physical Chemistry B, 104(15):3712–3720, 2000.
  * [49] A. Onufriev, D. A. Case, and D. Bashford.  Effective Born radii in the generalized Born approximation: the importance of being perfect.  Journal of Computational Chemistry, 23(14):1297–1304, 2002.
  * [50] A. Rahman.  Correlations in the Motion of Atoms in Liquid Argon.  Phys. Rev., 136:A405–A411, 1964.
  * [51] B. Roux and T. Simonson.  Implicit Solvent Models.  Biophys. Chem., 28:155–179, 1999.
  * [52] C. Sagui and T. A. Darden.  MOLECULAR DYNAMICS SIMULATIONS OF BIOMOLECULES: Long- Range Electrostatic Effects.  Annu. Rev. Biophys. Biomol. Struct., 28:155–179, 1999.
  * [53] M. F. Sanner, A. J. Olson, and J.-C. Spehner.  Reduced surface: an efficient way to compute molecular surfaces.  Biopolymers, 38(3):305–320, 1996.
  * [54] A. Scheraga, H. A.; Khalili, M.; Liwo.  Protein-Folding Dynamics: Overview of Molecular Simulation Techniques.  Annu. Rev. Phys. Chem., 58:57–83, 2007.
  * [55] T. Sutmann, G.; Gibbon, P.; Lippert.  Fast Methods for Long-Range Interactions in Complex Systems.  Forschungszentrum Jülich, 2011.
  * [56] H. Tjong and H. X. Zhou.  GBr6NL: A generalized Born method for accurately reproducing solvation energy of the nonlinear Poisson-Boltzmann equation.  Journal of Chemical Physics, 126:195102, 2007.
  * [57] J. Tomasi, B. Mennucci, and R. Cammi.  Quantum mechanical continuum solvation models.  Chem. Rev., 105:2999–3093, 2005.
  * [58] C. Wang and Y. Zhang.  Improving scoring-docking-screening powers of protein–ligand scoring functions using random forest.  Journal of computational chemistry, 38(3):169–177, 2017.
  * [59] E. Wang, H. Sun, J. Wang, Z. Wang, H. Liu, J. Z. Zhang, and T. Hou.  End-point binding free energy calculation with mm/pbsa and mm/gbsa: strategies and applications in drug design.  Chemical reviews, 119(16):9478–9508, 2019.
  * [60] J. Wang, H. Cao, J. Z. Zhang, and Y. Qi.  Computational protein design with deep learning neural networks.  Scientific reports, 8(1):1–9, 2018.
  * [61] J. Wang, C. H. Tan, Y. H. Tan, Q. Lu, and R. Luo.  Poisson-Boltzmann solvents in molecular dynamics Simulations.  Communications in Computational Physics, 3(5):1010–1031, 2008.
  * [62] Z. Wang, Y. Han, J. Li, and X. He.  Combining the fragmentation approach and neural network potential energy surfaces of fragments for accurate calculation of protein energy.  The Journal of Physical Chemistry B, 124(15):3027–3035, 2020.
  * [63] G.-W. Wei.  Protein structure prediction beyond alphafold.  Nature Machine Intelligence, 1:336–337, 2019.
  * [64] K. Wu and G. W. Wei.  Quantitative Toxicity Prediction Using Topology Based Multitask Deep Neural Networks.  Journal of Chemical Information and Modeling, 58:520–531, 2018\.
  * [65] K. Wu, Z. Zhao, R. Wang, and G. W. Wei.  TopP-S: Persistent Homology-Based Multi-Task Deep Neural Networks for Simultaneous Predictions of Partition Coefficient and Aqueous Solubility .  Journal of Computational Chemistry, 39:1444–1454, 2018.
  * [66] Y. C. Zhou, M. Feig, and G. W. Wei.  Highly accurate biomolecular electrostatics in continuum dielectric environments.  Journal of Computational Chemistry, 29(1):87–97, 2008.