# On the coupling instability of a gas jet impinging on a liquid film

David Barreiro-Villaverde1,2<EMAIL_ADDRESS>Anne Gosset3  Marcos
Lema3  Miguel A. Mendez2  1Universidade da Coruña, CITIC Research, Campus de
Elviña, 15071 A Coruña, Spain 2von Karman Institute for Fluid Dynamics, B-1640
Waterloosesteenweg 72, Sint-Genesius-Rode, Belgium 3Universidade da Coruña,
Campus Industrial de Ferrol, CITENI, 15403 Ferrol, Spain

###### Abstract

We investigate the dynamics of a gas jet impinging on a thin liquid film. This
configuration is relevant to the jet-wiping process and is unstable. In
particular, we complement previous works that focused on the wiping of liquids
with low Kapitza numbers (highly viscous liquids) by numerically analyzing the
wiping of liquids with much higher Kapitza numbers, more relevant to
industrial processes. The simulations are carried out by combining Volume of
Fluid (VOF) and Large Eddy Simulation (LES), and the dynamics of the gas-
liquid interaction is analyzed using extended multiscale Proper Orthogonal
Decomposition (emPOD). The resolution and flow details captured by the
simulations are unprecedented. The results show that, despite the vastly
different wiping conditions, the dynamics of the gas-liquid interaction is
remarkably similar. This opens new avenues to the study and the scaling of the
jet-wiping process.

###### keywords:

Jet wiping process, impinging gas jets, thin films

##  1 Introduction

Gas jets impinging on thin films are used in wiping processes in the coating
industry (Buchlin, 1997; Gosset & Buchlin, 2007). In these processes, the gas
jet acts as an air-knife which controls the thickness of the liquid film
deposited on the substrate. The uniformity of the final coating is known to be
limited by the instability of the gas-liquid interaction. Long-wave
disturbances were first revealed by Gosset (2007), numerically reproduced by
Myrillas et al. (2013), and extensively characterized by Gosset et al. (2019)
and Mendez et al. (2019b). These works have shown that the primary instability
of the wiping consists of large oscillations of the impinging jet, combined
with large waves in the liquid film, propagating both downstream and upstream
of the wiping region. An extensive analysis of this mechanism was presented by
Barreiro-Villaverde et al. (2021), who combined Volume of Fluid (VOF) and
Large Eddy Simulation (LES) with an extension of the multiscale Proper
Orthogonal Decomposition (mPOD, Mendez et al. 2019a) to study the interaction
between the two flows. This investigation revealed that the wiping instability
is essentially bidimensional, with waves originating as a result of both
displacement of the wiping region and modulation of its strength.

All the aforementioned studies focused on a narrow range of wiping conditions,
using low Kapitza number liquids:
$\mathrm{Ka}=\sigma\rho_{l}^{-1}\nu_{l}^{-4/3}g^{-1/3}$ (with $\rho_{l}$,
$\nu_{l}$ the liquid’s density and kinematic viscosity, $\sigma$ the gas-
liquid surface tension and $g$ the gravitational acceleration). These liquids
are characterized by a high viscosity and a low surface tension. It implies
that the film thickness is relatively large, and thus more “intrusive” for the
impinging gas jet. For instance, Gosset et al. (2019) and Barreiro-Villaverde
et al. (2021) focused on the wiping of viscous liquids with $\mathrm{Ka}$ in
the range $3-5$, typical of continuous painting processes or paper coating.
Nevertheless, the much more popular case of wiping in hot dip galvanization is
characterized by $\mathrm{Ka}\sim\mathcal{O}(10^{4})$ and significantly higher
Reynolds numbers in both the liquid film and the gas jet. A detailed
investigation of these conditions, however, is still out of the reach of both
experimental and numerical fluid dynamics. Studies of these conditions are
limited to theoretical models Hocking et al. 2011; Mendez et al. 2020a;
Ivanova et al. 2023 or 2D simulations Pfeiler et al. 2017b, a. The challenges
in high-fidelity simulations of the problem are discussed by Aniszewski et al.
(2020).

Therefore, an important open question is the extent to which the results
obtained for highly viscous fluids are relevant for the much broader spectra
of wiping conditions encountered in industry. This short article is an attempt
to tackle this question. Specifically, we complement the work in Barreiro-
Villaverde et al. (2021) with the investigation of the wiping of a liquid with
much higher $\mathrm{Ka}$. Although far from achieving full dynamic similarity
with galvanizing lines, as discussed in the following section, the
investigated conditions cover a completely different wiping regime than what
was previously reported in literature. Therefore, the results give an insight
into the dynamics and the scaling of the jet wiping instability, and perhaps
enable an educated extrapolation.

##  2 Selected Test Cases and Scaling Laws

We consider planar jet wiping, using a slot nozzle of opening $d_{n}$ and
width $W\gg d_{n}$, placed at a distance $Z_{n}$ from a vertical strip moving
upwards at a velocity $U_{p}$. The strip is flat, and the problem is treated
as isothermal. The relevant gas and liquid properties are $\rho_{g},\nu_{g}$
and $\rho_{l},\nu_{l}$, respectively, and $\sigma$. The nozzle’s stagnation
chamber is at a gauge pressure $\Delta P_{N}$, and the liquid thickness
downstream the wiping is $h_{f}$. The scaling of this configuration is
discussed by Gosset et al. 2019 and Mendez et al. 2020a. In addition to
$\mathrm{Ka}$, the relevant dimensionless numbers are the Reynolds numbers
$\mathrm{Re}_{f}=U_{p}h_{f}/\nu_{l}$ for the liquid film and
$\mathrm{Re}_{j}=U_{j}d/\nu_{g}$ for the gas jet, the dimensionless standoff
distance $\hat{Z}_{n}=Z_{n}/d_{n}$, and three numbers. Two of these link the
wiping “strength” of the impinging jet with the liquid properties. The first
is the wiping number ${\Pi}_{g}=\Delta P_{N}d_{n}/(\rho_{l}gZ_{n}^{2})$ and
relates the maximum pressure gradient (responsible for most of the wiping
work) to the liquid density, while the second is the shear number
$\mathcal{T}_{g}=\Delta P_{N}d_{n}/(Z(\rho_{l}g\mu_{l}U_{p})^{1/2})$ and
relates the maximum shear stress to the liquid viscosity (see Gosset et al.
2019). The third is the dimensionless final thickness
$\hat{h}_{f}=h_{f}/h_{0}$, with $h_{0}=\sqrt{\nu_{l}U_{p}/g}$ the maximum
thickness that can be withdrawn in the absence of wiping (see Spiers et al.
1974), and indicates the intensity of the wiping. Finally, the number
$h_{0}/Z_{n}$ measure the “intrusiveness” of the liquid film on the gas jet
flow, so as to characterize its geometrical confinement.

This article analyzes four wiping conditions, for which Table 1 collects the
aforementioned parameters. Two of these (Cases 1 and 2) are taken from
Barreiro-Villaverde et al. (2021) and correspond to experiments in Mendez et
al. (2019b). They consider the wiping of dipropylene glycol (DG) with a jet of
air ($\rho_{g}=1.2$ kg/m3 and $\nu_{g}=1.48\cdot 10^{-5}$ m2/s) with
$Z_{n}=18.5$ mm, $d_{n}=1.3$ mm and $U_{p}=0.34$ m/s and discharge velocity
$U_{j}\approx 26$ m/s (Case 1, with $\Delta P_{N}=425$ Pa) and $U_{j}\approx
38$ m/s (Case 2, with $\Delta P_{N}=875$ Pa). The other two cases (Cases 3 and
4) were added for this work and consider the wiping of water (W) by an air jet
with $Z_{n}=10$ mm, $d_{n}=1$ mm and $U_{p}=1$ m/s and discharge velocity
$U_{j}\approx 42$ m/s (Case 4, with $\Delta P_{N}=1$ kPa) and $U_{j}\approx
50$ m/s (Case 3, with $\Delta P_{N}=1,5$ kPa).

It is instructive to compare these wiping conditions with those in a
galvanizing line. The table reports the relevant parameters for a moderate
wiping of molten zinc by means of an air jet with $Z_{n}=15$ mm, $d_{n}=1.2$
mm and $U_{p}=3$ m/s, and discharge velocity $U_{j}\approx 160$ m/s ($\Delta
P_{N}\approx 20$ kPa). These would produce a final thickness of about $\approx
25\mu$m. The newly investigated conditions with water are significantly
different from previous test cases and closer to galvanizing conditions, as
evidenced in Table 1.

Table 1: Dimensional and dimensionless wiping conditions for cases 1 and 2 with dipropylene glycol (DG), with $\rho_{l}=1023$ kg/m3, $\nu_{l}=7.33\cdot 10^{-5}$ m2/s, $\sigma_{l}=0.032$ N/m, Cases 3 and 4 with water (W), with $\rho_{l}=1000$ kg/m3, $\nu_{l}=1\cdot 10^{-6}$ m2/s, $\sigma_{l}=0.073$ N/m, and an example of galvanizing conditions with zinc (Z), with $\rho_{l}=6500$ kg/m3, $\nu_{l}=4.5\cdot 10^{-7}$ m2/s, $\sigma_{l}=0.78$ N/m.  Liquid | Case ID | $\mathbf{h_{f}(\mu m)}$ | $\mathbf{h_{0}/Z_{n}}$ | $\mathbf{h_{f}/h_{0}}$ | $\mathbf{Re_{f}}$ | $\mathbf{Ka}$ |  $\mathbf{U_{j}}$ (m/s) | $\mathbf{Re_{j}}$ | $\mathbf{\Pi_{g}}$ | $\mathbf{\mathcal{T}_{g}}$
---|---|---|---|---|---|---|---|---|---|---
DG | 1 | 435 | 0.085 | 0.27 | 2.1 | 4.8 | 26.6 | 2280 | 0.16 | 0.24
2 | 263 | 0.16 | 1.4 | 38.2 | 3274 | 0.33 | 0.41
W | 3 | 27.2 | 0.032 | 0.085 | 27.2 | 3401 | 40.8 | 2690 | 1.02 | 2.23
4 | 21.3 | 0.067 | 21.3 | 50.0 | 3297 | 1.53 | 3.02
Z | 25 | 0.024 | 0.069 | 172 | 16440 | 160 | 12916 | 1.25 | 3.12

##  3 Methodology

The selected test cases are simulated in OpenFOAM, combining the algebraic VOF
formulation of the interFoam solver with the Smagorinski model for the
turbulence treatment of the gas jet flow. The numerical approach was
extensively validated for test cases 1 and 2 in Barreiro-Villaverde et al.
(2021), and is now applied to simulate a much more challenging set-up in cases
3 and 4. The domain and snapshots of the mesh (consisting of about 14 million
cells) are shown in figure 1. The domain includes the coating bath from which
the flat substrate is withdrawn, and a portion of the slot nozzle from which
the gas jet is released. It spans 50 mm above and below the jet axis in the
stream-wise direction ($x$), 20 mm in addition to the standoff distance
$Z_{n}$ in the cross-stream direction ($y$), and 20 mm in the span-wise
direction ($z$). The boundary conditions are indicated in figure 1, and the
gas jet is established through an inlet condition with prescribed stagnation
pressure $\Delta P_{N}$. The lateral patches (not shown) are set to cyclic.

Figure 1: Numerical domain, boundary conditions and mesh discretization.

The discretisation is particularly challenging because of the multiscale
nature of the problem. The cell size in the wiping region ($-Z_{n}<x<Z_{n}$)
is $\Delta_{x}=50\,\mu$m and $\Delta_{y}=2\,\mu$m across the film thickness.
This provides approximately 10-15 cells within the thickness of the final
film, and a minimum of 150 cells per dominant wavelength in the stream-wise
direction. The cell size in the span-wise direction is $\Delta_{z}=250\,\mu$m.
On the gas side, the LES index (defined as the ratio between the resolved and
the total turbulent kinetic energy) is kept above 80 % in the wiping region.
The simulations are initiated as in Barreiro-Villaverde et al. (2021) and
require $\approx$ 100 ms to move past the initial transient, after which 350
ms of physical time are simulated. The time step is $10^{-7}$ s to keep the
CFL number below 0.9, and a total of $n_{t}$=3500 snapshots are exported with
a sampling rate of 10 kHz. It requires approximately 1000 hours running in
parallel on 512 Intel E5–2680v3 CPUs from the Centro de Supercomputacion de
Galicia (CESGA) for one test case with water (case 3).

The results are processed using the extended mPOD as in Barreiro-Villaverde et
al. (2021). In particular, we decompose the film thickness $h(x,y,t)$ into
mPOD modes and identify the leading ones as those with the largest and paired
amplitudes (denoted as $\sigma_{r}$), having spatial ($\phi_{r}(x,y)$) and
temporal structures ($\psi_{r}(t)$) in quadrature. These correspond to
travelling waves. The mPOD modes are optimal modes for a prescribed frequency
partition (see Mendez et al. (2019a); Mendez (2023)) and have a band-limited
frequency content. The temporal structures of the leading waves in the liquid
film are then used to identify the most correlated coherent structures in the
gas jet. The latter are computed by projecting the velocity fields on the
temporal structures $\psi_{r}(t)$. The reader is referred to Barreiro-
Villaverde et al. (2021) for more details on the extended mPOD.

##  4 Results

We begin by examining the wave patterns in the liquid film in section 4.1, and
move to the study of the correlated jet structures in section 4.2. Finally,
section 4.3 discusses the interaction between the two phases. In all sections,
the new results on cases 3 and 4 are analysed along with the ones previously
obtained for cases 1 and 2.

###  4.1 Wave patterns on the liquid film

We describe the dynamics of the liquid film upstream the wiping region with
the help of figure 2, which compares case 1 (a) and case 3 (b). In each
figure, an instantaneous film thickness distribution is shown on the left,
with the liquid in red where $u>0$ and in green where $u<0$. On the right, it
is complemented with a contour plot of the stream-wise velocity component in
the midplane ($u(x,y,z=L_{z}/2)$), together with three velocity profiles taken
at the crest of a run-back wave. These profiles are (1) the one obtained from
the CFD computations, (2) the one obtained using the 2D Integral Boundary
Layer (IBL) model (equation 3.6 in Mendez et al. 2020a) with the flow rate,
film thickness and interface shear stress computed from CFD, and (3) the one
obtained with 2D IBL model neglecting the gas shear stress. An additional
movie of the film thickness distribution with a 2D velocity field taken at the
z-midplane is provided for cases 1 and 3.

In case 1, the liquid film is clearly bidimensional and features large waves
originating approximately $2.5$ mm below the jet axis at a frequency of about
$20$ Hz, corresponding to
$\hat{f}=f\mathrm{Ca}^{-1/3}(U_{p}\nu_{l}/g)^{-1/2}\approx 0.11$, with
$\mathrm{Ca}=\mu_{l}U_{p}/\sigma$ the capillary number. The reference
quantities are taken from the Shkadov-like scaling proposed in Gosset et al.
2019 for liquid films dragged by moving substrates. These falling waves have a
wavelength of the order of $\lambda\approx 12$ mm and evolve over an average
film thickness of the order of $h_{r}\approx 2h_{0}=$ 3 mm. The Reynolds
number in this region, defined as $Re_{r}=h_{r}\Delta U/\nu$, with $\Delta
U=\max(U_{p}-u|_{y=h})\approx 2U_{p}$ is of the order of $Re_{r}\approx 18$.
The flow is fairly laminar, and the bi-dimensional nature of the falling waves
is in line with what happens at the onset of interface stability in falling
liquid films. The flow pattern in the film displays a stagnation line at the
wavefront, but the assumption of a self-similar parabolic profile seems
appropriate. The contribution of the interface shear stress on the velocity
profile is much smaller than the viscous and gravitational one, as one might
expect from the low shear stress number in this case.

(a) Case 1.

(b) Case 3.

Figure 2:  Snapshots of the film dynamics in cases 1 (a) and 3 (b). On the
left, a 3D view of the liquid film colored in red where $u>0$ and in green
where $u<0$; and on the right, a velocity contour and streamline plot in the
midplane. The figure on the top compares the velocity profile at the wave
crest with the prediction of the IBL model with and without interface shear
stress. A movie of the 3D reconstruction of the liquid film and 2D velocity
field taken at the z-midplane is also provided for both cases.

Case 3 is clearly in a different wiping regime. The Reynolds number in this
region is of the order of $Re_{r}\approx 800$, and the waves rapidly undergo a
3D transition, at approximately $0.5$ mm below the jet axis. The velocity
field within the liquid is more complex, but several remarkable similarities
can be observed. First, the stream-wise velocity component maintains a
parabolic shape, even if the influence of the shear stress is more pronounced
than in case 1. Second, the frequency of wave formation ($f\approx$ 100-150
$Hz$) in the introduced scaling is $\hat{f}\approx$ 0.15-0.2, a value notably
close to the one in case 1. Third, and perhaps more importantly for the
following discussion, the waves are initially bi-dimensional in both cases.

(a) Case 1 (DG): $\hat{Z}=14.2$, $\Pi_{g}=0.16$, $\mathcal{T}_{g}=0.24$.

(b) Case 2 (DG): $\hat{Z}=14.2$, $\Pi_{g}=0.33$, $\mathcal{T}_{g}=0.41$.

(c) Case 3 (W): $\hat{Z}=10$, $\Pi_{g}=1.02$, $\mathcal{T}_{g}=2.23$.

(d) Case 4 (W): $\hat{Z}=10$, $\Pi_{g}=1.53$, $\mathcal{T}_{g}=3.02$.

Figure 3: Comparison between snapshots of the normalized film thickness
$\check{h}(\hat{x},\hat{z})$ (Original) with the leading wave patterns
detected via mPOD (mPOD filtered). The spectra of the temporal structure of
these modes are shown on the right as a function of the dimensional $f$ and
dimensionless $\hat{f}$ wave frequencies.

We now move to the liquid film downstream the wiping region. Figure 3 shows
snapshots of the liquid film thickness (left), its reconstruction using the
leading mPOD modes linked to the dominant travelling wave pattern (middle),
and the spectra of the associated temporal structures (right). A snapshot of
the four cases is included, and the wiping conditions are recalled in the
caption. Both the mPOD analysis and the plots are carried out using the
normalized film thickness, defined as
$\check{h}(x,z,t)=({h(x,z,t)-\overline{h}(x,z)})/{\sigma_{h}(x,z)}$ with
$\overline{h}(x,z)$ the average thickness and $\sigma_{h}(x,z)$ the thickness
standard deviation Barreiro-Villaverde et al. 2021. This normalization allows
the decomposition to equalize the importance of waves downstream and upstream
the wiping region despite the largely different thickness. The spectra
computed from the temporal structure of the leading modes are shown as a
function of the dimensional $f$ and dimensionless $\hat{f}$ wave frequencies.

In all the cases, the leading wave pattern is remarkably two-dimensional in
spite of the vastly different wiping conditions, and regardless of the three-
dimensional pattern arising upstream the wiping region with water. The
dominant wavelength upstream and downstream wiping differs because of the
different advection velocities in the two regions, but like in the run-back
film, the range of dimensionless frequency is surprisingly similar. In line
with the experimental findings on low $\mathrm{Ka}$ liquids in Gosset et al.
(2019), the frequency tends to increase with the wiping strength. On the other
hand, the gas jet Strouhal number $St=fZ_{n}/U_{j}$ based on the dominant
frequency of the waves, is in the range (0.01-0.04), well below the typical
values in free-jet instability and hydrodynamic feedback mechanisms Mendez et
al. 2018. This suggests a coupling between the gas jet and the liquid, and the
goal of the section that follows is to detect the gas jet structures linked to
those dominant wave patterns.

###  4.2 Gas jet structures correlated with the leading waves in the liquid
film

We now focus on the spatial structures in the gas jet flow using the extended
mPOD (emPOD). This technique allows revealing the flow structures that are
most correlated with the leading wave patterns in the liquid film (described
in figure 3). Figure 4 shows two representative snapshots of the gas velocity
field $\bm{u}=(u,v)$ and its projection for case 1 (a) and case 3 (b). Cases 2
and 4 lead to similar results and are thus omitted for brevity, but an
animation of the film thickness evolution and gas jet flow for each test case
is provided in the supplementary material. In each figure, the contour on the
left is the mean centred flow field (i.e.
$\bm{u}^{\prime}=\bm{u}(\bm{x})-\overline{\bm{u}}$) while the one on the right
is the projection of the flow field on the leading mPOD modes of the
normalized film thickness.

The original flow fields reveal the complexity of the gas jet dynamics. After
impinging the liquid interface, the jet flow splits into two “side jets”
evolving along the liquid interface. The velocity gradient on the sides of the
gas jet core triggers the formation of vortices (label “$v_{s}$”) that
destabilize the jet near the impingement area and induce small-scale
disturbances at the liquid film interface (label “$v_{f}$”), which are rapidly
damped due to viscuous and capillary damping Barreiro-Villaverde et al. 2023.
This dynamics, however, is not correlated with the wavy patterns and thus, is
not visible in the projected fields on the right. Here, the leading structures
consist of a larger vortex induced by the falling liquid waves (label
“$v_{r}$”) and a rigid oscillation of the impinging jet (label “$d$”). In the
case of low Kapitza liquids, this mechanism was extensively documented in
Myrillas et al. (2013); Mendez et al. (2019b); Barreiro-Villaverde et al.
(2021), where a hypothesis on their interaction with the liquid film dynamics
was formulated. In particular, it was postulated that the jet oscillation was
linked to the vortex-liquid film interaction. The current results with high
$\mathrm{Ka}$ liquids show that this interaction is much weaker in these
conditions because of the comparatively smaller thickness of the liquid film.
Yet, it is shown here for the first time that these structures “$v_{r}$” and
“$d$” persist, suggesting that the underlying mechanism of the unstable
interaction between the gas jet and the liquid film is the same.

(a) Case 1 (DG)

(b) Case 3 (W)

Figure 4:  Mean-centred $u^{\prime}$ and emPOD projected fields for cases 1
(a) and 3 (b). The plots are scaled such that the horizontal axis ranges
between $0<y/Z<1$, and the vertical axis between $-10<x/d<10$. A movie for
each test case in table 1 is also provided.

###  4.3 Dynamics of the gas-liquid instability

Finally, the gas jet-liquid film interaction is analyzed by means of the
spatio-temporal contours of the dimensionless pressure gradient
$\partial_{\hat{x}}\hat{p}=\partial_{x}p/(\rho_{l}g)$ (top) and normalized
film thickness $\check{h}$ (bottom) for cases 1 (left) and 3 (right) in figure
5. Both quantities are taken at the mid-plane. The stream-wise coordinate is
normalized using the jet opening $d_{n}$ while the time scale is scaled as in
Mendez et al. (2020b), i.e. $\hat{t}=tCa^{1/3}(U_{p}\nu_{l}/g)^{1/2}$. The
evolution of the pressure gradient reveals the amplitude of the wiping
unsteadiness. The latter is significantly stronger in case 1, resulting in a
much larger modulation of the pressure gradient and, thus, of the wiping
strength. This difference is clearly linked to the different levels of
“intrusiveness” of the wiped liquid.

The dynamics of the liquid film is more regular in case 1 than in case 3.
Upstream the wiping region, the characteristic lines of the waves are gently
curved, but continuous in case 1, while in case 3, waves occasionally coalesce
or reverse direction. The space-time contours at the bottom are complemented
with a plot of the temporal evolution of (1) the impingement point (dashed red
line) and (2) the wiping point (continuous white line), i.e. the location of
the highest pressure gradient in $x<0$. In all cases, it appears that the
amplitude of the jet oscillation has a minor role compared to the large
fluctuations experienced by the wiping point. These fluctuations are
characterized by a “slow” downward dynamic and a “fast” upward dynamic, as if
the mechanisms governing these two stages were completely different. The
downward shift of the wiping point occurs at the scale of the wave advection
time in the falling liquid waves, while the upward shift occurs at the scale
of the gas jet advection time. This corroborates the hypothesis that this
mechanism is linked to an interaction between the gas jet and the liquid film.
The most remarkable result of this analysis is that the same dynamics are
observed in both high $\mathrm{Ka}$ and low $\mathrm{Ka}$ cases, even if the
wiping conditions and intrusiveness of the film are radically different.
Because the dimensionless frequencies of the waves remain within the same
range, we may infer that the interaction is mostly driven by the liquid film.
Considering that the dimensionless groups related to the liquid film dynamics
are in fair similarity between water and zinc (table 1), these results suggest
that coupling dynamics observed in water might not differ significantly from
the one occurring in galvanization. It is striking that the wavy defects
observed on galvanized products have typical wavelengths in the range 10-15
mm, which corresponds to dimensionless frequencies of $\hat{f}=(0.11-0.165)$,
totally in line with what is found here for largely different conditions.

Case 1 (DG)

Case 3 (W)

Figure 5: Spatio-temporal contours of the dimensionless pressure gradient
(top) and normalized film thickness (bottom) for cases 1 (left) and 3 (right).
The thickness contours are complemented with the temporal evolution of the
impingement point (dashed red line) and the wiping point (continuous white
line), i.e. the location of the highest pressure gradient in $x<0$.

##  5 Conclusions

We have numerically investigated the two-phase coupling instability taking
place between an impinging gas jet and a liquid film dragged by an upward
moving substrate. The two-phase CFD computations cover largely different
wiping conditions with a level of detail that is unprecedented in literature.

It is found that the waves emerge two-dimensional in the impingement region in
all the investigated conditions. In some cases, the dominant 2D patterns
undergo a 3D transition due to intrinsic instabilities in the run-back flow,
or due to the impact of small-scale vortices in the final film. From the gas
jet side, the correlation analysis reveals two main structures acting at the
time scale of the leading wave pattern: a symmetric oscillation around the jet
axis (pattern “$d$”), and a deflection of the lower side jet triggered by the
periodic formation of waves (pattern “$v_{r}$”). It is shown that the second
mechanism has a stronger impact on the pressure gradient and, thus, on the
formation of the waves.

Finally, in spite of the very different flow regimes analyzed in this work, it
is remarkable that the dynamics of the gas-liquid interaction is qualitatively
similar, and that the wave frequency scales reasonably well using a purely
liquid-based scaling. Although the system locks at a certain frequency that
depends on both the gas jet and liquid film, these results suggest a dominant
role of the liquid film in the coupling.

## References

  * Aniszewski et al. (2020) Aniszewski, W., Saade, Y., Zaleski, S. & Popinet, S. 2020 Planar jet stripping of liquid coatings: Numerical studies. International Journal of Multiphase Flow 132, 103399.
  * Barreiro-Villaverde et al. (2023) Barreiro-Villaverde, D., Gosset, A., Lema, M. & Mendez, M. A. 2023 Damping of three-dimensional waves on coating films dragged by moving substrates. Physics of Fluids 35 (7), 72110.
  * Barreiro-Villaverde et al. (2021) Barreiro-Villaverde, D., Gosset, A. & Mendez, M.A. 2021 On the dynamics of jet wiping: Numerical simulations and modal analysis. Physics of Fluids 33 (6).
  * Buchlin (1997) Buchlin, J. M. 1997 Modelling of gas jet wiping in Thin Liquid Films and Coating Processes. In VKI Lecture Series. Rhode-Saint-Genese.
  * Gosset (2007) Gosset, A. 2007 Study of the interaction between a gas flow and a liquid film entrained by a moving surface. PhD thesis, Université Libre de Bruxelles.
  * Gosset & Buchlin (2007) Gosset, A. & Buchlin, J. M. 2007 Jet Wiping in Hot-Dip Galvanization. Journal of Fluids Engineering 129 (4), 466.
  * Gosset et al. (2019) Gosset, A., Mendez, M. A. & Buchlin, J. M. 2019 An experimental analysis of the stability of the jet wiping process: Part I – Characterization of the coating uniformity. Experimental Thermal and Fluid Science 103, 51–65.
  * Hocking et al. (2011) Hocking, G. C., Sweatman, W. L., Fitt, A. D. & Breward, C. 2011 Deformations during jet-stripping in the galvanizing process. Journal of Engineering Mathematics 70 (1-3), 297–306.
  * Ivanova et al. (2023) Ivanova, T., Pino, F., Scheid, B. & Mendez, M. A. 2023 Evolution of waves in liquid films on moving substrates. Physics of Fluids 35 (1), 013609.
  * Mendez (2023) Mendez, M A 2023 Generalized and Multiscale Modal Analysis, pp. 153–181. Cambridge University Press.
  * Mendez et al. (2019a) Mendez, M. A., Balabane, M. & Buchlin, J. M. 2019a Multi-scale proper orthogonal decomposition of complex fluid flows. Journal of Fluid Mechanics 870, 988–1036.
  * Mendez et al. (2019b) Mendez, M. A., Gosset, A. & Buchlin, J.-M. 2019b Experimental analysis of the stability of the jet wiping process, part II: Multiscale modal analysis of the gas jet-liquid film interaction. Experimental Thermal and Fluid Science 106, 48–67.
  * Mendez et al. (2020a) Mendez, M. A., Gosset, A., Scheid, B., Balabane, M. & Buchlin, J. M. 2020a Dynamics of the jet wiping process via integral models, arXiv: 2004.13400.
  * Mendez et al. (2020b) Mendez, M. A., Hess, D., Watz, B. B. & Buchlin, J. M. 2020b Multiscale proper orthogonal decomposition (mPOD) of TR-PIV data-a case study on stationary and transient cylinder wake flows. Measurement Science and Technology 31 (9).
  * Mendez et al. (2018) Mendez, M. A., Scelzo, M. T. & Buchlin, J. M. 2018 Multiscale modal analysis of an oscillating impinging gas jet. Experimental Thermal and Fluid Science 91, 256–276.
  * Myrillas et al. (2013) Myrillas, K., Rambaud, P., Mataigne, J. M., Anderhuber, M., Gardin, P., Vincent, S. & Buchlin, J. M. 2013 Numerical modeling of gas-jet wiping process. Chemical Engineering and Processing: Process Intensification 68, 26–31.
  * Pfeiler et al. (2017a) Pfeiler, C., Eßl, W., Reiss, G., Ecker, W., Riener, C. K. & Angeli, G. 2017a LES-VOF Simulation and POD Analysis of the Gas-Jet Wiping Process in Continuous Galvanizing Lines. Steel Research International 88 (9), 1600507.
  * Pfeiler et al. (2017b) Pfeiler, C., Eßl, W., Reiss, G., Riener, C. K, Angeli, G. & Kharicha, A. 2017b Investigation of the gas-jet wiping process- two-phase large eddy simulations elucidate impingement dynamics and wave formation on zinc coatings. Steel Research International 88 (9), 1600507.
  * Spiers et al. (1974) Spiers, R. P., Subbaraman, C. V. & Wilkinson, W. L. 1974 Free coating of a Newtonian liquid onto a vertical surface. Chemical Engineering Science 29 (2), 389–396.


# On the importance of catalyst-adsorbate 3D interactions for relaxed energy
predictions

Alvaro Carbonero
Mila
<EMAIL_ADDRESS>
Alexandre Duval
Mila, CentraleSupélec
&Victor Schmidt
Mila, Université de Montréal
&Santiago Miret
Intel labs
&Alex Hernandez-Garcia
Mila, Université de Montréal &Yoshua Bengio
Mila, Université de Montréal &David Rolnick
Mila, McGill University

###### Abstract

The use of machine learning for material property prediction and discovery has
traditionally centered on graph neural networks that incorporate the geometric
configuration of all atoms. However, in practice not all this information may
be readily available, e.g. when evaluating the potentially unknown binding of
adsorbates to catalyst. In this paper, we investigate whether it is possible
to predict a system’s relaxed energy in the OC20 dataset while ignoring the
relative position of the adsorbate with respect to the electro-catalyst. We
consider SchNet, DimeNet++ and FAENet as base architectures and measure the
impact of four modifications on model performance: removing edges in the input
graph, pooling independent representations, not sharing the backbone weights
and using an attention mechanism to propagate non-geometric relative
information. We find that while removing binding site information impairs
accuracy as expected, modified models are able to predict relaxed energies
with remarkably decent MAE. Our work suggests future research directions in
accelerated materials discovery where information on reactant configurations
can be reduced or altogether omitted.

##  1 Introduction

Materials discovery is a strong driver of innovation, unlocking new materials
with tailored properties that benefit various domains including energy
efficiency, transportation systems and electronics (Butler et al., 2018). Yet,
it faces significant hurdles. Characterizing new materials is computationally
expensive (Oganov et al., 2019), even when replacing lab experiments by
quantum mechanical simulations like the Density Functional Theory (DFT) (Kohn
et al., 1996). Besides, exploration of materials is hindered by the vastness
of the search space, which encompasses myriad compositions, atomic
arrangements and properties (Pyzer-Knapp et al., 2022).

To overcome these challenges, researchers have turned to Machine Learning (ML)
for accelerated materials discovery (Zhang et al., 2023; Miret et al., 2022;
Lee et al., 2023; Song et al., 2023) for two primary reasons: First, ML holds
the power to quickly model materials’ properties (i.e., to evaluate
candidates) using Geometric Graph Neural Networks (GNNs). Second, generative
ML can automatically propose new and consistent material candidates.

In recent years, ML has emerged as a crucial tool for the discovery of
electro-catalysts, which play a key role in promoting renewable energy
processes and sustainable chemical production, including the production of
ammonia for fertilizers and hydrogen (Zitnick et al., 2020). The Open Catalyst
Project (OCP) (Chanussot et al., 2021) has significantly contributed to this
field by releasing an extensive dataset of 1,281,040 DFT relaxations of
catalyst-adsorbate pairs, selected from a pool of meanigful candidates111The
adsorbate refers to the molecule involved in the electrochemical reaction that
is accelerated through the introduction of a catalyst material, represented by
a semi-infinite periodic substructure.. This dataset was specifically designed
to train ML models to predict the relaxed energy of 3D adsorbate-catalyst
(adslab) systems, a critical property that influences a catalyst’s activity
and selectivity (Tran and Ulissi, 2018), or its effectiveness for a specific
chemical reaction222In other words, it quantifies the extent to which the
catalyst reduces the energy required for a chemical reaction to occur.. OCP
has facilitated significant advancements in catalysis discovery, with ML
models increasingly bridging the performance gap with DFT simulations, while
offering a speed advantage of several orders of magnitude (Musaelian et al.,
2022; Gasteiger et al., 2022; Passaro and Zitnick, 2023).

While research has primarily focused on predicting material properties (Wieder
et al., 2020), it’s equally crucial to efficiently explore the vast search
space of potential catalyst candidates (Wei et al., 2019; Zhang et al., 2022).
This involves generating consistent candidates automatically, a challenge due
to the intricate process of creating adslab samples. This process, further
detailed in Appendix A, includes positioning an adsorbate molecule (e.g., H2O)
with a catalyst in 3D space, meaning cutting a surface through the catalyst
bulk, selecting the adsorbate’s spatial orientation, and sampling a plausible
binding site (Chanussot et al., 2021). These steps, which are time-consuming
and challenging to model, determine the input configuration of the adslab
sample and significantly influence the relaxed energy prediction (Lan et al.,
2022). An additional challenge is that this process and its relevance to the
actual in-lab efficiency of a real material is not yet fully understood by
chemists (Deshpande et al., 2020).

In light of these challenges, this paper explores the possibility of
predicting the relaxed energy of an adslab without co-locating the adsorbate
and the catalyst in the same 3D space. This direction of study is valuable for
several reasons: (1) to better understand the role of adsorbate-catalyst
geometry in determining the relaxed energy; (2) to reduce reliance on the
exact input configuration, which often correlates with local energy minima,
(3) to avoid the complexity and computational cost of determining a good input
configuration (e.g., finding a binding site and adsorbate orientation). Our
primary objective is quantify the loss in accuracy that occurs when the
geometric relationship between inputs is unavailable.

To achieve this, we propose four modifications of existing GNN architectures,
collectively referred to as Disconnected GNNs. All four modifications enable
the base architecture to make relaxed energy predictions without explicitly
modeling geometric interactions between the adsorbate and the catalyst. We
then evaluate the trade-off of omitting these interactions through experiments
on the OC20 dataset and suggest future directions to overcome the limits
associated with the input adslab configuration.

##  2 Methods

Our methods test the assumption that meaningful results can be obtained in
predicting the relaxed energy even when the geometric interactions between the
adsorbate and the catalyst are unknown. Specifically, we propose four
Disconnected GNN models. These models can leverage as their backbone any
underlying GNN model that predict properties of 3D atomic systems using the
general pipeline detailed under “Refresher” below, where we make targeted
changes to the pipeline. In this paper, we investigate the effect of the
proposed Disconnected GNN architecture on three backbone models: SchNet
(Schütt et al., 2017), DimeNet++ (Klicpera et al., 2020) and FAENet (Duval et
al., 2023).

Refresher. Recall that GNNs applied to 3D atomic systems mostly use the
following pipeline. (1) Graph creation: construct a graph representation of
the input point cloud systems (e.g. adslab), represented with atomic numbers
and 3D atom positions. (2) Embedding: derive node/edge embeddings based on
atomic numbers and geometric information (e.g. relative atom positions). (3)
Interaction blocks: apply a fixed number of message passing layers (Gilmer et
al., 2017) to update the node and/or the edge embeddings, using geometric
information and preserving data symmetries. (4) Output block: project final
atom representations into scalar values (e.g. atom-properties). For graph-
level prediction, a global pooling (e.g. sum) is performed to aggregate atom
predictions.

(a) Baseline.

(b) Disconnected baseline.

(c) Independent pooling.

(d) Independent Backbones.

(e) Attention.

Figure 1: Architectures are presented in order of increasing complexity. 1(a)
illustrates the standard property prediction pipeline for Geometric GNNs, with
no modifications. In 1(b), edges between parts are removed. In 1(c), the
pooling step is performed separately on atom representations for the adsorbate
and the catalyst. Resulting embeddings are concatenated and passed through an
MLP. In 1(d), each part has its own model with distinct weights. In 1(e), an
interaction block with artificial edges and weights allows node embeddings
between both parts to communicate.

###  2.1 Disconnected baseline

The baseline Disconnected GNN model does not create any edges between the
catalyst and the adsorbate in the graph creation step. By removing all such
edges, no relative geometric information will pass from the adsorbate to the
catalyst when modelling the system, since GNNs only propagate information
through graph edges. Since this modification only relates to the input data
representation, this change can be applied to any backbone GNN architecture,
like FAENet, SchNet or DimeNet++ for instance. It is illustrated in Figure
1(b).

Note that we still use the original OC20 dataset where the adsorbate/catalyst
atom positions have been determined after selecting a binding site and an
adsorbate orientation. But since GNNs use atom relative positions and not
absolute atom positions to preserve translation equivariance, disconnecting
the adsorbate-catalyst graph is enough to discard totally the input
configuration. Processing such a disconnected graph is equivalent to
considering both components independently. Associated results thus inform
about the importance of the relative adslab geometric information to each
model.

###  2.2 Independent pooling

This method builds upon the disconnected baseline defined above, also removing
all edges between the adsorbate and the catalyst. As illustrated in Figure
1(c), there are three additional changes. First, the projection block now
outputs a hidden representation $\mathbf{h}\in\mathbb{R}^{\left\lfloor
H/2\right\rfloor}$ for each node instead of a scalar quantity, where $H$ is
the dimension of the embedding. Second, the resulting node representations are
pooled separately for each component, leading to distinct hidden
representations for the adsorbate and the catalyst. Lastly, we concatenate
these two representations and pass them to an additional dense layer to
compute the energy of the system. This final MLP layer gives more expressive
power to the model than the disconnected baseline because we explicitly model
non-geometric interactions between the adsorbate and the catalyst.

###  2.3 Independent backbones

This approach, illustrated in Figure 1(d), builds upon independent pooling.
The main difference is that catalysts and adsorbates are assigned distinct GNN
models (hence “independent”) instead of sharing backbone weights. The
motivation is that the adsorbate and the catalyst have very different sizes
and roles. Hence, it might be beneficial to have independent GNNs, where each
one is trains for its particular input. Similarly to independent pooling, both
models are modified to produce a hidden representation of their respective
parts. They are then concatenated and passed through an MLP.

###  2.4 Attention

This model (see Figure 1(e)) builds upon the independent model by allowing the
nodes of both graphs to communicate during the interaction layers.

We modify the graph creation process to produce a heterogeneous graph where we
add the following weighted edges to previous models’ edges. For every catalyst
node $i$ and adsorbate node $j$, we create the new edges $(i,j)$ and $(j,i)$
with weights $z$ and $-z$ respectively, where $z$ is the $z$-axis coordinate
of $x_{i}$. Besides, in the embedding step, we pass the edge weights through
an RBF layer and a linear layer to increase their dimension. Notice that this
model can run without locating the adsorbate and catalyst in the same 3D plane
since no information regarding the adsorbate’s location (e.g. distance,
relative position) is given. The motivation behind these weighted edges is to
make the model aware of the closeness of a node to the catalyst’s surface,
while respecting the restrictions of the setup. As noted in (Duval et al.,
2022), fixed atoms not near the surface of the catalyst (tag 0 nodes) can be
removed without hurting efficiency, indicating that closeness to surface plays
a factor.

We now describe how this graph is used. Let $F$ denote a “normal” intra-
interaction layer, let $h^{l}_{\text{ads}}$ and $h^{l}_{\text{cat}}$ denote
the node embeddings of the adsorbate and the catalyst at interaction layer
$l$, and let $e_{\text{ads}}$ and $e_{\text{cat}}$ stand for the edges of the
adsorbate and catalyst, respectively. To produce $h^{l+1}_{\text{ads}}$ and
$h^{l+1}_{\text{cat}}$, we first calculate
$h^{\prime}_{\text{ads}}=F(h^{l}_{\text{ads}},e_{\text{ads}})$ and
$h^{\prime}_{\text{cat}}=F(h_{\text{cat}},e_{\text{cat}})$. Subsequently, we
calculate
$h^{\prime}=\text{GATConv}(\text{Concat}(h^{\prime}_{\text{ads}},h^{\prime}_{\text{cat}}),e^{\prime})$
where $e^{\prime}$ is the set of new weighted edges and GATConv is a
convolution layer as described in (Veličković et al., 2017). Finally,

$h^{l+1}_{\text{ads}},h^{l+1}_{\text{cat}}=\text{Norm}(h^{\prime}_{\text{ads}}+h^{l}_{\text{ads}}),\text{Norm}(h^{\prime}_{\text{cat}}+h^{l}_{\text{cat}}).$

##  3 Evaluation

Dataset: We use the OC20 IS2RE (Zitnick et al., 2020) dataset, which involves
the direct prediction of the relaxed adsorption energy from the initial atomic
structure, i.e. a graph regression task requiring E(3)-invariance. OC20 comes
with $450K$ training samples and a predefined train/val/test split. We
evaluate models on the four distinct splits of the validation set ($\sim 25K$
samples each): In Domain (ID), Out of Domain Adsorbates (OOD-ads), Out of
Domain catalysts (OOD-cat), and Out of Domain Adsorbates and catalysts (OOD-
both).

Metrics: We measure accuracy on each validation split through the energy MAE.
Running time is measured with the throughput at inference time, i.e. the
average number of samples per second that a model can process in its forward
pass, on similar GPU types.

Baselines: We compare base results of the original models (considering 3D
adsorbate-catalyst interactions) with the four disconnected architectures
detailed in Section 2. Note that we apply PhAST components (Duval et al.,
2022) on all models to improve performance and inference time.

Figure 2: Comparison of MAE performance drops across our four proposed
modifications to three base architectures (left), and their associated impact
on model throughput (right). We can see the the disconnected baseline
consistently achieves the best performance. In terms of throughput, it is on
par with other modifications for DimeNet++ and FAENet, while it is slower than
all other modifications for SchNet (except for the Attention modification).

Results. From Figure 2, we see that the independent model architecture
performs best among all proposed disconnected GNNs. However, it presents a
performance drop in MAE of 20%, 31%, and 28% with respect to baseline
architectures, for Schnet, DimeNet++ and FAENet respectively. A complete
numerical table of results is reported in Table B.

Performance is significantly worse for disconnected GNNs than for GNNs with
access to the adslab configuration; however, these results are promising in
several ways.

Firstly, a drop in performance was expected since we neglect important
information by omitting to modelise the adsorbate-catalayst 3D interaction.
The important take-away is that our disconnected GNN manages to produce decent
energy predictions while being agnostic to the adslab’s input configuration.
Future work will probably manage to reduce further this performance gap. Note
that we discuss some axes of improvements in Appendix C.

Secondly, the OC20 dataset often contains multiple input configuration samples
for the same catalyst-adsorbate pair, each with a different relaxed energy
target. This cannot be captured by our disconnected GNN, leading to multiple
targets for the same sample, which can impair training dynamics. Furthermore,
when evaluating on a dataset that does not present this feature, e.g. the 10k
split of OC20 IS2RE, disconnected GNNs match (or outperform) baseline models,
which is very promising. We present and briefly discuss these results in
Appendix D.

The OC20-Dense dataset (Lan et al., 2022) was recently released, containing
multiple input configurations (e.g. binding site, adsorbate orientation) for
each adsorbate-catalyst pair along with the associated relaxed energy.
Building upon this dataset, we will be able to train disconnected GNN models
to predict the minimum relaxed energy of each adslab over all possible input
configurations. The obvious benefit of such approach is that it marginalises
over adslab binding sites and adsorbate orientations, avoiding the need to
explore all configurations to actually find the relaxed energy global minima
of this adsorbate-catalyst pair, which is expected to be the main reaction
driver in real life experiments (Lan et al., 2022).

##  4 Conclusion

In this paper, we have explored predicting an adsorbate-catalyst system’s
relaxed energy without co-locating them in 3D space. We proposed four
disconnected models built on top of existing Geometric GNNs architectures, and
showed through experiments that this task was possible with a manageable
performance loss. Looking forward, we see the potential of using such
disconnected models to find the global energy minima of adslab systems while
circumventing the need to consider all possible input adslab configurations.
This sounds particularly promising for generative methods and should
contribute to the acceleration of catalysis discovery.

## References

  * Ashcroft and Mermin (2022) Neil W Ashcroft and N David Mermin.  _Solid state physics_.  Cengage Learning, 2022.
  * Boes et al. (2019) Jacob R Boes, Osman Mamun, Kirsten Winther, and Thomas Bligaard.  Graph theory approach to high-throughput surface adsorption structure generation.  _The Journal of Physical Chemistry A_ , 123(11):2281–2285, 2019.
  * Butler et al. (2018) Keith T Butler, Daniel W Davies, Hugh Cartwright, Olexandr Isayev, and Aron Walsh.  Machine learning for molecular and materials science.  _Nature_ , 559(7715):547–555, 2018.
  * Chanussot et al. (2021) Lowik Chanussot, Abhishek Das, Siddharth Goyal, Thibaut Lavril, Muhammed Shuaibi, Morgane Riviere, Kevin Tran, Javier Heras-Domingo, Caleb Ho, Weihua Hu, et al.  Open catalyst 2020 (oc20) dataset and community challenges.  _ACS Catalysis_ , 11(10):6059–6072, 2021.
  * Deshpande et al. (2020) Siddharth Deshpande, Tristan Maxson, and Jeffrey Greeley.  Graph theory approach to determine configurations of multidentate and high coverage adsorbates for heterogeneous catalysis.  _npj Computational Materials_ , 6(1):79, 2020\.
  * Duval et al. (2022) Alexandre Duval, Victor Schmidt, Santiago Miret, Yoshua Bengio, Alex Hernández-García, and David Rolnick.  Phast: Physics-aware, scalable, and task-specific gnns for accelerated catalyst design.  _arXiv preprint arXiv:2211.12020_ , 2022.
  * Duval et al. (2023) Alexandre Duval, V. Schmidt, A. Garcia, Santiago Miret, Fragkiskos D. Malliaros, Y. Bengio, and D. Rolnick.  Faenet: Frame averaging equivariant gnn for materials modeling.  _International Conference on Machine Learning_ , 2023.  doi: 10.48550/arXiv.2305.05577.
  * Gasteiger et al. (2022) Johannes Gasteiger, Muhammed Shuaibi, Anuroop Sriram, Stephan Günnemann, Zachary Ulissi, C Lawrence Zitnick, and Abhishek Das.  Gemnet-oc: developing graph neural networks for large and diverse molecular simulation datasets.  _arXiv preprint arXiv:2204.02782_ , 2022.
  * Gilmer et al. (2017) J. Gilmer, S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural message passing for quantum chemistry.  _International Conference On Machine Learning_ , 2017.
  * Jain et al. (2013) Anubhav Jain, Shyue Ping Ong, Geoffroy Hautier, Wei Chen, William Davidson Richards, Stephen Dacek, Shreyas Cholia, Dan Gunter, David Skinner, Gerbrand Ceder, et al.  The materials project: A materials genome approach to accelerating materials innovation, apl mater.  2013\.
  * Klicpera et al. (2020) Johannes Klicpera, Shankari Giri, Johannes T Margraf, and Stephan Günnemann.  Fast and uncertainty-aware directional message passing for non-equilibrium molecules.  _Preprint arXiv:2011.14115_ , 2020.
  * Kohn et al. (1996) Walter Kohn, Axel D Becke, and Robert G Parr.  Density functional theory of electronic structure.  _The Journal of Physical Chemistry_ , 100(31):12974–12980, 1996.
  * Kresse and Hafner (1994) Georg Kresse and Jürgen Hafner.  Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium.  _Physical Review B_ , 49(20):14251, 1994.
  * Lan et al. (2022) Janice Lan, Aini Palizhati, Muhammed Shuaibi, Brandon M. Wood, Brook Wander, Abhishek Das, Matt Uyttendaele, C. Lawrence Zitnick, and Zachary W. Ulissi.  Adsorbml: Accelerating adsorption energy calculations with machine learning.  _arXiv preprint arXiv: Arxiv-2211.16486_ , 2022.
  * Lee et al. (2023) Kin Long Kelvin Lee, Carmelo Gonzales, Marcel Nassar, Matthew Spellings, Mikhail Galkin, and Santiago Miret.  Matsciml: A broad, multi-task benchmark for solid-state materials modeling, 2023.
  * Miret et al. (2022) Santiago Miret, Kin Long Kelvin Lee, Carmelo Gonzales, Marcel Nassar, and Matthew Spellings.  The open MatSci ML toolkit: A flexible framework for machine learning in materials science.  _Preprint arXiv:2210.17484_ , 2022.
  * Musaelian et al. (2022) Albert Musaelian, Simon Batzner, Anders Johansson, Lixin Sun, Cameron J Owen, Mordechai Kornbluth, and Boris Kozinsky.  Learning local equivariant representations for large-scale atomistic dynamics.  _Preprint arXiv:2204.05249_ , 2022.
  * Oganov et al. (2019) Artem R Oganov, Chris J Pickard, Qiang Zhu, and Richard J Needs.  Structure prediction drives materials discovery.  _Nature Reviews Materials_ , 4(5):331–348, 2019\.
  * Ong et al. (2013) Shyue Ping Ong, William Davidson Richards, Anubhav Jain, Geoffroy Hautier, Michael Kocher, Shreyas Cholia, Dan Gunter, Vincent L Chevrier, Kristin A Persson, and Gerbrand Ceder.  Python materials genomics (pymatgen): A robust, open-source python library for materials analysis.  _Computational Materials Science_ , 68:314–319, 2013.
  * Passaro and Zitnick (2023) Saro Passaro and C Lawrence Zitnick.  Reducing so (3) convolutions to so (2) for efficient equivariant gnns.  _arXiv preprint arXiv:2302.03655_ , 2023.
  * Pyzer-Knapp et al. (2022) Edward O Pyzer-Knapp, Jed W Pitera, Peter WJ Staar, Seiji Takeda, Teodoro Laino, Daniel P Sanders, James Sexton, John R Smith, and Alessandro Curioni.  Accelerating materials discovery using artificial intelligence, high performance computing and robotics.  _npj Computational Materials_ , 8(1):84, 2022\.
  * Schütt et al. (2017) Kristof Schütt, Pieter-Jan Kindermans, Huziel Enoc Sauceda Felix, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  _Advances in neural information processing systems_ , 30, 2017.
  * Song et al. (2023) Yu Song, Santiago Miret, and Bang Liu.  MatSci-NLP: Evaluating scientific language models on materials science language tasks using text-to-schema modeling.  In _Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 3621–3639, Toronto, Canada, July 2023. Association for Computational Linguistics.  doi: 10.18653/v1/2023.acl-long.201.  URL https://aclanthology.org/2023.acl-long.201.
  * Tran and Ulissi (2018) Kevin Tran and Zachary W Ulissi.  Active learning across intermetallics to guide discovery of electrocatalysts for co2 reduction and h2 evolution.  _Nature Catalysis_ , 1(9):696–703, 2018.
  * Veličković et al. (2017) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  _arXiv preprint arXiv:1710.10903_ , 2017.
  * Wei et al. (2019) Jing Wei, Xuan Chu, Xiang-Yu Sun, Kun Xu, Hui-Xiong Deng, Jigen Chen, Zhongming Wei, and Ming Lei.  Machine learning in materials science.  _InfoMat_ , 1(3):338–358, 2019.
  * Wieder et al. (2020) Oliver Wieder, Stefan Kohlbacher, Mélaine Kuenemann, Arthur Garon, Pierre Ducrot, Thomas Seidel, and Thierry Langer.  A compact review of molecular property prediction with graph neural networks.  _Drug Discovery Today: Technologies_ , 37:1–12, 2020.
  * Zhang et al. (2023) Mengchun Zhang, Maryam Qamar, Taegoo Kang, Yuna Jung, Chenshuang Zhang, Sung-Ho Bae, and Chaoning Zhang.  A survey on graph diffusion models: Generative ai in science for molecule, protein and material.  _arXiv preprint arXiv:2304.01565_ , 2023.
  * Zhang et al. (2022) Xu Zhang, Yun Tian, Letian Chen, Xu Hu, and Zhen Zhou.  Machine learning: a new paradigm in computational electrocatalysis.  _The Journal of Physical Chemistry Letters_ , 13(34):7920–7930, 2022.
  * Zitnick et al. (2020) C Lawrence Zitnick, Lowik Chanussot, Abhishek Das, Siddharth Goyal, Javier Heras-Domingo, Caleb Ho, Weihua Hu, Thibaut Lavril, Aini Palizhati, Morgane Riviere, et al.  An introduction to electrocatalyst design using machine learning for renewable energy storage.  _Preprint arXiv:2010.09435_ , 2020.

##  Appendix A Adslab generation

This section summarises the adslab creation process of the OC20 dataset. Refer
to the original paper [Chanussot et al., 2021] for more details.

Adsorbate and catalyst surface selection. The first step is to select the
adsorbate and the catalyst that will compose the adslab. For the adsorbate, it
extremely simple: it is sampled randomly from a set of 82 molecules that are
chosen for their utility to renewable energy applications. For the catalyst
surface, the process is divided in three stages. First, we choose the number
of distinct chemical elements composing it. It can be a unary material (5%
chance), a binary material (65%) or a ternary material (30%). These elements
are chosen from a set of 55 elements comprising reactive nonmetals, alkali
metals, transition metals, etc. Next, a stable bulk material is randomly
selected from the 11K samples of the Materials Project [Jain et al., 2013]
with the number of elements chosen in the first step. Lastly, all
symmetrically distinct surfaces333a surface is a cut through the bulk operated
using Miller indices from the material with Miller indices [Ashcroft and
Mermin, 2022] less than or equal to 2 are enumerated, including possibilities
for different absolute positions of surface plane. From this list of surfaces,
one is randomly selected.

Input adslab configuration. The objective of the second step is to place the
adsorbate and the catalyst surface in the same 3D plane. Using pymatgen’s
Voronoi tesselation algorithm [Ong et al., 2013] on surface atoms and
adsorbate possible binding sites, Catkit [Boes et al., 2019] enumerates a list
of symmetrically distinct adsorption sites along with suggested per-site
adsorbate orientations. From this list, an adsorption configuration is
randomly selected, yielding the initial adslab structure. This initial
structure is then relaxed using Density Functional Theory simulations,
performed using the Vienna Ab Initio simulation Package (VASP) [Kresse and
Hafner, 1994].

##  Appendix B Full results table

Baseline / MAE | ID | OOD-ad | OOD-cat | OOD-both | Average | Throughput (samples/s)
---|---|---|---|---|---|---
5pt.  SchNet | $0.631$ | $0.687$ | $0.625$ | $0.626$ | $\boldsymbol{0.642}$ | $16001\pm 2264$
Base disconnected-SchNet | $0.940$ | $1.037$ | $0.932$ | $0.947$ | $0.964$ | $18863\pm 3034$
Independent pooling-SchNet | $0.766$ | $0.847$ | $0.762$ | $0.789$ | $0.791$ | $18590\pm 2230$
Independent models-SchNet | $0.761$ | $0.805$ | $0.764$ | $0.764$ | $\underline{0.773}$ | $12220\pm 1593$
Attention-SchNet | $0.779$ | $0.852$ | $0.762$ | $0.797$ | $0.797$ | $5652\pm 876$
DimeNet++ | $0.577$ | $0.693$ | $0.568$ | $0.621$ | 0.615 | $631\pm 125$
Base disconnected-DimeNet++ | $0.922$ | $1.059$ | $0.914$ | $0.989$ | $0.971$ | $667\pm 143$
Independent pooling-DimeNet++ | $0.781$ | $0.0.992$ | $0.767$ | $0.901$ | $0.860$ | $650\pm 139$
Independent models-DimeNet++ | $0.767$ | $0.884$ | $0.754$ | $0.818$ | $\underline{0.806}$ | $374\pm 68$
Attention-DimeNet++ | $0.777$ | $0.952$ | $0.763$ | $0.870$ | $0.841$ | $329\pm 57$
5pt.  FAENet | $0.554$ | $0.623$ | $0.546$ | $0.578$ | 0.575 | $1263\pm 766$
Base disconnected-FAENet | $0.884$ | $0.937$ | $0.878$ | $0.870$ | $0.892$ | $1549\pm 1027$
Independent pooling-FAENet | $0.698$ | $0.813$ | $0.697$ | $0.766$ | $0.743$ | $1553\pm 1049$
Independent models-FAENet | $0.690$ | $0.782$ | $0.689$ | $0.760$ | $\underline{0.731}$ | $803\pm 327$
Attention-FAENet | $0.732$ | $0.814$ | $0.715$ | $0.756$ | $0.754$ | $902\pm 565$
5pt. |  |  |  |  |  |

Table 1: MAE and inference time for various GNNs and their disconnected
counterpart on OC20 IS2RE, averaged over 3 runs. Average MAE is computed over
all validation splits. Best results are shown in bold, second best underlined.
Overall, disconnected models show a lower accuracy than baselines but they are
able to predict relaxed energy relatively well despite omitting geometric
interaction between parts in the adslab.

##  Appendix C Axes of improvement

There are two main ways in which Disconnected GNN models can be further
improved. The first one is to improve upon dense layers as a way to combine
the final hidden representations of the adsorbate and the catalyst. The second
one is to find a way to effectively allow the nodes of the adsorbate and
catalyst to communicate, and to do this while ommitting the relative position
and orientation between the adsorbate and the catalyst. This is attempted in
the attention models through the inter-interaction layers. Given that they do
not consistently outperform the independent backbone models, it is likely that
the way in which attention models use GAT convolution layers hinders the
model’s prediction. We still hope that this can be fixed since the attention
models can outperform baseline models when trained on the 10k dataset, as we
described in Appendix D.

##  Appendix D Results in the 10k dataset

Only 89% of the OC20 IS2RE dataset consists of adsorbate-catalyst systems
which are unique up to the combination of adsorbate id, bulk id, and bulk
cell. That is, to Disconnected GNNs, 11% of the dataset consists of inputs
with multiple targets due to distinct binding sites, which impairs training
dynamics. When compared to the OC20 10k IS2RE dataset, every single adsorbate-
catalyst system has a unique combination of adsorbate id, bulk id, and bulk
cell. Thus, evaluating disconnected GNNs on the 10k dataset is a useful way to
assess their performance on a suitable data tailored to them. This is
corroborated by the results in Table D.

Baseline 10k/ MAE | ID | OOD-ad | OOD-cat | OOD-both | Average | Throughput (samples/s)
---|---|---|---|---|---|---
5pt.  SchNet | $1.151$ | $1.245$ | $1.104$ | $1.213$ | $1.179$ | $18602\pm 4373$
Base disconnected-SchNet | $1.246$ | $1.236$ | $1.214$ | $1.144$ | $1.210$ | $25487\pm 3693$
Independent pooling-SchNet | $0.962$ | $1.038$ | $0.935$ | $0.961$ | $0.974$ | $24832\pm 3041$
Independent models-SchNet | $0.976$ | $0.985$ | $0.946$ | $0.935$ | $\boldsymbol{0.961}$ | $\underline{16501\pm 3283}$
Attention-SchNet | $0.979$ | $1.072$ | $0.961$ | $0.968$ | $0.995$ | $\boldsymbol{7952\pm 1261}$
DimeNet++ | $0.851$ | $0.969$ | $0.802$ | $0.875$ | 0.874 | $623\pm 123$
Base disconnected-DimeNet++ | $1.044$ | $1.078$ | $1.026$ | $0.971$ | $1.023$ | $649\pm 140$
Independent pooling-DimeNet++ | $0.908$ | $1.022$ | $0.888$ | $0.928$ | $0.937$ | $752\pm 114$
Independent models-DimeNet++ | $0.889$ | $0.971$ | $0.876$ | $0.883$ | $\underline{0.905}$ | $\underline{366\pm 71}$
Attention-DimeNet++ | $0.889$ | $1.017$ | $0.865$ | $0.915$ | $0.921$ | $\boldsymbol{330\pm 60}$
5pt.  FAENet | $1.003$ | $1.017$ | $0.992$ | $1.004$ | $1.004$ | $\underline{1417\pm 792}$
Base disconnected-FAENet | $1.104$ | $1.144$ | $1.118$ | $1.061$ | $1.107$ | $1492\pm 1085$
Independent pooling-FAENet | $1.022$ | $1.009$ | $0.999$ | $0.938$ | $0.992$ | $1626\pm 1123$
Independent models-FAENet | $0.891$ | $1.020$ | $0.871$ | $1.040$ | $\underline{0.955}$ | $1089\pm 737$
Attention-FAENet | $0.894$ | $0.924$ | $0.863$ | $0.872$ | $0.888$ | $\boldsymbol{916\pm 567}$
5pt. |  |  |  |  |  |

Table 2: MAE and inference time for various GNNs and their disconnected
counterpart on OC20 10k IS2RE, averaged over 3 runs. Average MAE is computed
over all validation splits. Disconnected models all show lower accuracy. In
every model, every subsurface catalyst atom (tag 0 atoms) is removed as
suggested in [Duval et al., 2022] to improve inference time.

From Table D, we see that independent models perform on par or better (for
FAENet) than baseline models when trained on the 10k dataset. In particular,
we see a performance improvement of 18%, -4%, and 5% for SchNet, DimeNet++,
and FAENet respectively. Different than in the all dataset, Attention-FAENet
outperforms Independent models-FAENet by 7.5%. These results are a sign that
Disconnected GNNs can perform well when tested against the data for which they
were designed for. The OC20-Dense dataset [Lan et al., 2022] is thus an
exciting development as it will allow for the creation of the datasets that
these models need.