A scalable and quantum-accurate foundation model for biomolecular force fields via linearly tensorized quadrangle attention

Domain: physics
Published: 2026-03-07
Editorial outcome: minor_revision

Abstract
Accurate atomistic biomolecular simulations are vital for understanding disease mechanisms and
drug discovery, yet existing methods struggle to balance quantum-mechanical accuracy with
computational scalability. Classical force fields often lack precision, while quantum methods are
computationally prohibitive for complex biological systems. Here we show that LiTEN, a scalable
equivariant neural network, resolves this dilemma by efficiently modeling complex three- and four-
body interactions with linear complexity via Linearly Tensorized Quadrangle Attention. We
introduce LiTEN-FF, a foundation model pre-trained on extensive datasets to ensure broad chemical
generalization across diverse molecular spaces. We demonstrate that LiTEN achieves state-of-the-art
accuracy on standard benchmarks, consistently outperforming leading approaches in both precision
and speed. Furthermore, LiTEN-FF enables comprehensive modeling tasks, ranging from geometry
optimization to free energy surface construction, with high computational efficiency for large
biomolecules. This framework provides a physically grounded, versatile foundation for advanced
biomolecular modeling and drug design applications.

Introduction
Accurate atomistic simulations are essential for understanding disease mechanisms, advancing drug
discovery, and designing biomaterials by elucidating biomolecular structures and functions1, 2.
Traditional methods fall into two categories: classical force fields and quantum chemical calculations.
Classical force fields enable large-scale molecular dynamics (MD) simulations due to their high
computational efficiency but struggle to accurately model bond rearrangements, transition state
energetics, and fine conformational nuances due to fixed-parameters , often leading to discrepancies
in dynamic simulations3. In contrast, quantum chemical methods explicitly describe electronic structures and bond transformations, making them invaluable for
reaction pathway analysis and small molecule optimization4, 5. However, their high computational
cost restricts their use to short-timescale or small-system simulations, limiting their applicability to
complex biomolecular dynamics. Thus, achieving a practical balance between physical accuracy and
computational efficiency remains a fundamental and ongoing challenge in biomolecular modeling6,
7.
Recent advances in deep learning have opened new venues for addressing challenges in
biomolecular modeling. Machine learning interatomic potentials (MLIPs), trained on quantum
chemical data, have demonstrated a favorable trade-off between accuracy and computational
efficiency, enabling physically grounded modeling of complex molecular systems8, 9, 10. The
evolution of MLIPs follows three key trends: (1) transitioning from global molecular descriptors to
localized atomic environment representations; (2) advancing from approximate symmetry handling
to rigorous physical equivariances; and (3) extending from simple two-body to comprehensive many-
body interactions11, 12. This progression has given rise to several influential methodological
frameworks. SchNet13 pioneered continuous-filter convolutional networks for end-to-end prediction
of molecular energies and forces. PaiNN14 incorporated rotation-equivariant message passing,
significantly enhancing the model’s capacity to encode molecular symmetries. NequIP15 and
MACE16 developed graph neural network architectures that strictly comply with E(3) equivariance,
employing high-order tensor products to model many-body interactions, achieving state-of-the-art
(SOTA) performance on small molecules. AIMNet217 proposed hierarchical multi-scale interaction
modules to capture long-range electronic polarization effects with robust generalizability. The ANI18
series focused on drug-like molecules, combining local atomic environment descriptors with large-
scale training datasets for practical drug screening workflows.
Despite the promise of MLIPs in various applications, critical challenges remain in biomolecular
modeling. First, while datasets have expanded in scale and coverage, they are fragmented and rely
on ab initio methods with varying accuracy (e.g., DFT19, MP220, CCSD(T)21), lacking a unified
computational benchmark22, 23. This methodological heterogeneity introduces biases in potential
energy surfaces (PES), hindering cross-dataset training and transfer learning. Second, while datasets
such as QM924, ANI-1x, and AIMNet2 cover the general organic chemical space, their representation
of key biomolecular subspaces (e.g., carbohydrates, nucleotides, metal complexes, and polar residues)
remains insufficient, limiting model generalizability and predictive stability in real biological systems.
Given the limitations in existing datasets, current models often struggle with limited
generalization and reduced prediction accuracy. To improve performance, many approaches have
adopted more complex neural network architectures, which significantly increase computational
inference cost25. High-order equivariant models (e.g., NequIP and MACE) demonstrate superior
physical consistency and many-body modeling but suffer from low inference efficiency due to
complex construction and propagation of high-order tensors, especially for medium-to-large systems.
In contrast, lower-order models like AIMNet2, built upon Cartesian coordinates without explicit
tensor products, achieve faster inference speed, making them more suitable for industrial
deployment. However, they sacrifice predictive accuracy for non-local coupling or large systems, as
their inability to represent high-order tensors restricts capturing subtle many-body effects critical for
local potential energy surface curvature and coupling directions. Therefore, in the design of MLIPs,
completely omitting tensor product operations within the Cartesian coordinate framework can
substantially boost computational efficiency, though this comes at the cost of enforcing symmetry
constraints and compromising high-order many-body modeling. Conversely, employing high-order
tensor products enhances representational power but incurs significant computational costs,
limiting practicality for large molecular systems. Consequently, there is an urgent need for a
framework that balances physical consistency, modeling accuracy, and computational efficiency,
alongside new solutions for pretraining data and strategies.
In this work, we present LiTEN, an equivariant tensor message-passing model built in Cartesian
coordinates that combines physical consistency, expressive power, and inference efficiency. Its core
innovation is the Tensorized Quadrangle Attention (TQA) mechanism, which unifies the modeling
of three-body and four-body interactions centered on edges and nodes with linear complexity. This
enables efficient representation of many-body couplings in complex molecular conformations.
Unlike traditional equivariant networks that rely on explicit high-order tensor products and
spherical harmonic expansions, LiTEN reparameterizes high-order tensor structures using dot and
cross products between vectors, compactly encoding interactions while preserving equivariance. This
approach significantly reduces computational graph complexity and circumvents the performance
bottleneck associated with Clebsch-Gordan coefficient calculations. This design not only improves
the model’s structural awareness in dense many-body fields but also offers a physically motivated and
scalable solution for AI-driven force field modeling in bioorganic systems. Evaluations on multiple
benchmark datasets, including rMD1726, MD2227 and Chignolin28, demonstrate that LiTEN
achieves SOTA performance on nearly all test subsets, outperforming established methods such as
MACE, NequIP, and Allegro29. Built on the LiTEN architecture, we introduce LiTEN-FF, a robust
foundation model for atomistic simulation. To balance quantum-level accuracy and broad chemical
generalization, LiTEN-FF is trained through a two-stage strategy: pretraining on the large-scale
nablaDFT30 dataset (~16 million conformations of drug-like molecules with elements H, C, N, O,
S, Cl, F, Br), and then fine-tuned on the higher-accuracy SPICE31 dataset (~2 million bio-organic
molecules covering a wider element range including metals and halogens). This hierarchical
approach equips LiTEN-FF with strong transferability across diverse chemical and biological systems,
including solvated environments. In terms of modeling capability, LiTEN-FF achieves quantum
accuracy in conformer optimization across systems of varying sizes, with an average RMSD of just
0.048 Å compared to DFT references. In MD, it reproduces bond length and angle distributions
with KL divergences as low as 0.001–0.01. On the TorsionNet20632 benchmark, it outperforms
existing MLIPs and medium-level DFT methods, achieving a MAE of only 0.19 kcal/mol relative to
CCSD(T). Its predicted free energy surfaces closely match experimental observations, confirming its
suitability for thermodynamic modeling. In periodic aqueous simulations, LiTEN-FF accurately
reproduces radial distribution functions and performs similarly to MACE-OFF33 in modeling
peptide free energy landscapes. For computational performance, LiTEN-FF delivers up to a 10×
speedup over leading models like MACE-OFF on large systems (~1000 atoms), while maintaining
or exceeding their accuracy. It further supports batch-mode conformer search within MLIPs, scaling
up to 40× faster with increased batch size—highlighting its strength in high-throughput molecular
modeling. In summary, LiTEN and LiTEN-FF overcome key limitations of existing force fields,
delivering a scalable, quantum-accurate solution that pushes the boundaries of MLIPs in both speed
and accuracy. Their integration with structure prediction and molecular design pipelines promises
to accelerate real-world applications in drug discovery and materials science.

Results
Overview of LiTEN-FF
LiTEN adopts a hierarchical architecture (Fig. 1A), consisting of an embedding layer designed to
capture local atomic geometric features, followed by three core modules responsible for tensor-based
message passing, edge-based vector cross-product interactions, and node-based vector-scalar fusion.
The architecture culminates with output heads that predict total molecular energy and atomic forces.
Importantly, the atomic forces are computed as the negative gradient34 of the energy with respect to
atomic coordinates, ensuring physical consistency and stability for MD simulations.
In the model design, each atom 𝑖 is assigned a first-order equivariant vector 𝑢⃗⃗⃗ , Given the
𝑖
neighbor set 𝒩(𝑖) centered on atom 𝑖 , we perform a weighted aggregation of the normalized
displacement vectors 𝑢⃗⃗⃗⃗⃗ from all neighbors to obtain a first-order geometric descriptor:
𝑖𝑗
𝒖⃗⃗⃗ = ∑ 𝒖⃗⃗⃗⃗⃗ , 𝒖⃗⃗⃗⃗⃗ = 𝒓⃗⃗⃗⃗ 𝒊 ⃗ 𝒋 . (1)
𝒊 𝑗∈𝒩(𝑖) 𝒊𝒋 𝒊𝒋 |𝒓⃗⃗⃗⃗ 𝒊 ⃗ 𝒋|
This vector 𝑢⃗⃗⃗ encodes the dominant directional characteristics (i.e., first-order anisotropy) of
𝑖
the atom within its local environment. Based on the inner product of directional vectors, a three-
body structural embedding can be constructed without explicit angle calculations:
|𝒖⃗⃗⃗ |𝟐 = ∑ 𝒖⃗⃗⃗⃗⃗ ⋅𝒖⃗⃗⃗⃗⃗⃗ = ∑ cos𝜃 . (2)
𝒊 𝑗,𝑘∈𝒩(𝑖) 𝒊𝒋 𝒊𝒌 𝑗,𝑘 𝑗𝑖𝑘
This quadratic form is equivalent to the sum of the cosines of the angles between adjacent edges.
When neighboring atoms tend to align collinearly, this value reaches its maximum, reflecting the
angular consistency of the atom’s local environment. To better distinguish conformations with
different symmetries or principal axis distributions, we introduce a cubic composite tensor with two
main purposes. Firstly, it aims to embed angular information along with the magnitudes of atomic
vectors into the node scalars, thereby enhancing the model’s sensitivity to directional shifts in the
local conformation. Its expression is given by:
|𝒖⃗⃗⃗ |3 = (∑cos𝜃 )⋅|𝒖⃗⃗⃗ |. (3)
𝒊 I𝑗𝑖𝑘 𝒊
Secondly, this cubic term integrates both self-interaction and cross-interaction information,
effectively capturing weakly coupled many-body structures. Its approximate expansion is as follows:
|𝒖⃗⃗⃗⃗ |3 = (𝐴T+𝐵)3/2 ≈ 𝐴3/2+ 3 𝐴1/2𝐵+ 3 𝐴−1/2𝐵2. (4)
𝒊 2 8
Here, we define 𝐴 = ∑ |𝑢⃗⃗⃗⃗⃗ | 2 as the self-interaction term, and 𝐵 = ∑ 𝑢⃗⃗⃗⃗⃗ ⋅𝑢⃗⃗⃗⃗⃗ as the cross-
𝑗 𝑖𝑗 𝑗≠𝑘 𝑖𝑗 𝑖𝑘
interaction term. The above formula is a Taylor expansion valid under the condition |𝐵/𝐴| ≪ 1.
This approximation relies on the selection of neighboring atoms. Typically, local neighbors are
obtained within a cutoff radius, which comprehensively covers the atomic distribution around the
target atom. In the ideal case of a uniform atomic distribution, 𝐵 ≈ 0 , and in other cases, 𝐵
generally remains small, effectively describing weakly coupled many-body interactions35, 36. When the
local atomic arrangement exhibits a distinct directional bias (e.g., linear or L-shaped conformations),
this term’s response is enhanced, whereas it approaches zero in perfectly symmetric conformations
(e.g., tetrahedral), thereby improving the distinguishability of directional structures. These invariants,
constructed from directional tensors, collectively form a multibody descriptive basis that does not
require manual construction of angular terms.
To model higher-order conformational constraints, such as torsional potentials and spatial
rigidity, we propose an edge-oriented four-body interaction mechanism. For a quadruple (𝑖,𝑗,𝑘,𝑙),
where each pair of atoms is connected through shared edges, forming a chain-like structure, we
define the cosine of the dihedral (torsion) angle between adjacent edges as follows:
(𝒖⃗⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )⋅(𝒖⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )
𝒊𝒌 𝒊𝒋 𝒊𝒋 𝒋𝒍
cos𝜙 = . (5)
|𝒖⃗⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ |⋅|𝒖⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ |
𝒊𝒌 𝒊𝒋 𝒊𝒋 𝒋𝒍
This formula is equivalent to the cosine of the angle between the normal vectors of two adjacent
planes, capturing the spatial torsion information present in the four-atom chain.
By incorporating the sine of the angle, sin𝜃, and considering that the angular component cos𝜃
is already captured in the three-body interactions, we construct the following rotation-invariant four-
body geometric descriptor:
∑ sin𝜃 sin𝜃 cos𝜙 = ∑(𝒖⃗⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )⋅(𝒖⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )
𝑘𝑖𝑗 𝑙𝑗𝑖 I𝒊𝒌 𝒊𝒋 𝒊𝒋 𝒋𝒍
𝑘,𝑙∈𝒩(𝑖,𝑗) 𝑘,𝑙 (6)
= (𝒖⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )⋅(𝒖⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗⃗ ).
𝒊 𝒊𝒋 𝒊𝒋 𝒋
This information is encoded into edge features within the model and used for calculating
attention weights in subsequent iterations. Specifically, the weights in the implicit attention
mechanism37 are computed by directly summing the embeddings of nodes and edges, followed by
applying activation and scaling, and then combining with a distance decay factor:
⊤
𝛼 = [ ϕ( ℎ + ℎ + 𝑒 ) ⋅ 𝛼 ] ⋅R(𝑑 ). (7)
𝑖𝑗 𝑖 𝑗 𝑖𝑗 𝑖𝑗
This module significantly enhances the model’s sensitivity to molecular conformational
distortions, rigid fragments, and deviations from coplanarity, allowing it to more accurately capture
complex molecular structures such as conjugated systems, helical chains, and protein backbones. It
is worth noting that VisNet constructs dihedral angles using the vector rejection method. Although
this approach can reflect local geometric relationships to some extent, it does not explicitly normalize
the magnitude of the rejection vectors when aggregating neighboring vectors, which may weaken the
geometric interpretability of the resulting expressions. In contrast, our formulation based on cross-
product tensors provides a more direct correspondence to real geometric relationships and endows
the descriptors with clearer physical meaning. Furthermore, our model extends the four-body
interactions to higher-order many-body interactions and incorporates coupling relationships among
different interaction orders, thereby maintaining consistent physical representation in more
complex molecular systems.
In the LiTEN framework, all geometric and physical descriptors are constructed through tensor
operations within atomic neighborhoods, enabling the three-body and four-body modules to achieve
linear computational complexity O(N), which is significantly more efficient than traditional
combinatorial approaches. At the same time, we intentionally maintain simplicity in the design of
the message-passing modules, employing lightweight structures such as implicit attention
mechanisms to avoid unnecessary complexity. On this basis, LiTEN maintains excellent predictive
accuracy and numerical stability. Overall, LiTEN combines high computational efficiency with
strong generalization ability, ensuring robust applicability across diverse chemical environments
without relying on any predefined atom types. As shown in Fig. 1B-1C, the foundation model
LiTEN-FF, pretrained on NablaDFT and SPICE datasets, powers LiTEN-FF’s broad applications,
including vacuum-phase conformation optimization, free energy surface (FES) modeling, dihedral
scanning, and accurate geometric property prediction. Moreover, the model supports efficient
conformer search and is compatible with periodic systems. In summary, LiTEN-FF strikes an optimal
balance between accuracy, scalability, and generality, offering a powerful tool for molecular
simulations, AI-driven force fields, and large-scale molecular modeling.
Tests of model accuracy and generalization for quantitative properties
To systematically evaluate the predictive performance and generalization ability of the LiTEN model
across varying molecular scales, we selected two representative standard benchmark datasets: rMD17
and MD22. These datasets cover typical small and large molecular systems, respectively, offering a
comprehensive assessment of the model’s performance across diverse scenarios in terms of structural
complexity, atomic count, and chemical diversity.
The rMD17 dataset consists of multiple representative small organic molecules such as Aspirin,
Benzene, and Ethanol. These molecules generally have low molecular weights, relatively rigid
conformations, and stable electronic structures and geometries. They are commonly used to test the
model’s sensitivity and accuracy in capturing subtle atomic interactions within low-dimensional,
rigid systems. For these molecules, energy differences mainly arise from minor conformational
perturbations, demanding very high accuracy in force prediction. As a result, this dataset serves as
an important benchmark for evaluating the model’s fine-scale structural modeling capabilities at the
small molecular scales. In contrast, the MD22 dataset primarily contains more complex, biologically
relevant large molecular systems such as the tripeptide Ac-Ala3-NHMe, carbohydrates like Stachyose,
DNA fragments (AT-AT, AT-AT-CG-CG), and supramolecular structures such as the Buckyball
Catcher and Double-walled Nanotube (Dw_Nanotube). These systems feature high atom counts,
numerous degrees of freedom, and rich non-covalent interactions including hydrogen bonding, π-π
stacking, and van der Waals forces. MD22 effectively tests the model’s expressive power, scalability,
and structural generalization ability in high-dimensional structural spaces. Therefore, its provides a
stringent benchmark to determine whether LiTEN can robustly generalize to complex molecular
systems. 
For comparison, we selected several mainstream advanced DL-driven force field models, which
can be roughly divided into two groups. The first group is based on spherical harmonics convolution,
such as MACE, NequIP, and BOTNet38. These methods explicitly consider rotational equivariance
in tensor construction, making them highly suitable for capturing direction-sensitive angles and
many-body interactions. The second group relies on on Cartesian coordinate message passing
mechanisms, such as PaiNN, Allegro, and VisNet39, which maintain responsiveness to directional
information while generally offering higher efficiency and structural adaptability. Additionally, we
included several representative models, namely GemNet40, ACE41, ANI, and sGDML42 as references
to cover a broader range of paradigm designs and modeling strategies.
In the experimental setup, we uniformly used the optimal hyperparameter configurations
recommended in the literature for each model to ensure fairness in comparison. Except for VisNet
and MACE, which used their recommended numbers of layers (9 layers and 2 layers respectively),
other models employed a 6-layer network architecture with a hidden dimension size of 256. During
the training process, all models were run under the same training-validation-test splits, and standard
error metrics (MAE) were used to evaluate the energy (in meV or kcal/mol) and force (in meV/Å or
kcal/mol/Å) predictions.
On the rMD17 dataset, LiTEN achieved the best or second-best results in the vast majority of
small molecule systems, as shown in Table 1. Notably, for well-known molecules such as aspirin,
benzene, and ethanol, its force prediction errors were significantly lower than those of other models,
demonstrating outstanding performance in high-precision modeling of small perturbations. This
result also highlights the efficiency and stability of the model’s directional aggregation mechanism
and equivariant tensor modeling in capturing subtle conformational changes. On the MD22 dataset,
LiTEN showed excellent performance in several structurally complex large molecule systems, as
shown in Table 2. For example, in systems like Ac-Ala3-NHMe, Stachyose, AT-AT, and
Dw_Nanotube, LiTEN significantly outperformed most comparative models in terms of energy and
force prediction errors. Especially in scenarios with a large number of atoms, its prediction errors
exhibited only negligible fluctuations as the system size increased. This fully demonstrates that
LiTEN not only has good scalability when modeling large-scale molecules but also possesses
structural robustness and strong generalization ability across different molecular types.
Overall, as a many-body graph neural network with equivariant tensor modeling capabilities,
LiTEN achieves near-physical-limit high-precision predictions in small molecule systems and also
demonstrates excellent modeling capability and stability in large molecule systems. Compared to
existing state-of-the-art models, LiTEN shows significant advantages in terms of energy and force
prediction accuracy across multiple test sets and molecular structures. This validates the effectiveness
of its proposed directional aggregation, tensorized four-body attention, and implicit structure-aware
mechanisms, reflecting strong physical consistency and broad adaptability. In the future, LiTEN is
expected to play a pivotal role in complex tasks such as drug design, biomolecular modeling, and
new material prediction. It has the potential to emerge as one of the efficient and general-purpose
deep learning force field models.
Training Efficiency Evaluation on the Chignolin Dataset
Chignolin is a mini-protein composed of 10 amino acid residues, characterized by a well-defined β-
hairpin structure43. Due to its compact size and rich conformational diversity, it is widely used as a
benchmark dataset for evaluating protein folding, MD simulations, and the performance of machine
learning force field models. In this study, we adopted the Chignolin dataset introduced in the VisNet
work to evaluate both the training speed and prediction accuracy of various models, with a particular
focus on practical efficiency at moderate molecular scales44. All models were trained using a batch size
of 4.
The evaluation results demonstrate that the LiTEN family of models achieves significantly faster
training speeds while maintaining high prediction accuracy. Specifically, LiTEN-128 (with a hidden
dimension of 128) achieved an average training time of 22.75 ms per molecule, ranking just behind
PaiNN and ET in speed, while maintaining the lowest peak memory usage among all tested models
at 3912 MiB. The higher-capacity LiTEN-256 model, with a hidden dimension of 256, further
improved the prediction accuracy to 0.3641 kcal/mol/Å, with a training time of 31.21 ms per
molecule and memory consumption of 13949 MiB, still substantially lower than those of
mainstream models such as Allegro (33965 MiB) and NequIP (23841 MiB).
As shown in Fig. 2, the LiTEN series demonstrates a well-balanced trade-off among training
speed, memory utilization, and prediction accuracy. In particular, LiTEN-128 exhibits strong
potential for deployment in high-throughput simulations and resource-limited environments,
attributed to its low hardware demands and fast training. Overall, the results on the realistic protein
system Chignolin validate LiTEN’s advantages in modeling many-body interactions, highlighting its
robust generalization and practical applicability.
LiTEN-FF construction via hierarchical two-stage training on large-scale quantum
chemical datasets
To leverage the full potential of the wealthy large-scale quantum chemical data, we adopted a two-
stage hierarchical training strategy. In the initial stage, we performed large-scale pretraining on the
nablaDFT dataset, aiming to comprehensively capture the diverse electronic structure patterns
prevalent in drug-like molecules. For data partitioning, we reserved geometry optimization
trajectories as the test set and selected a subset of samples with significant scaffold diversity and
conformational perturbations for hyperparameter tuning. The remaining conformations were used
for model training (see Supplementary Table 1-2), laying a solid foundation for the model’s
generalization in downstream tasks. Specifically, the model architecture consists of six equivariant
convolutional layers, each with 256-dimensional hidden channels, and a cutoff radius of 5 Å, which
effectively covers the typical physical range of most intramolecular chemical interactions. The
model’s equivariance property ensures that the predicted energies and forces adhere to physical
symmetries under rotation and translation, thereby improving both its generalization ability and
physical consistency. The model was trained using energies in Hartree and forces in Hartree/Å, with
a weighted MAE loss. After 100 epochs of training, the model achieved excellent predictive accuracy
on the nablaDFT test set, with MAEs of 0.32 milliHatree for energy and 0.20 milliHatree/Å for
force, equivalent to approximately 0.20 kcal/mol and 0.126 kcal/mol/Å, respectively. These error
levels approach those of high-accuracy quantum chemical methods such as DFT, and notably fall
well within the widely accepted threshold (~1 kcal/mol) of what is known as chemical accuracy,
demonstrating the potential of LiTEN-FF Nab as a quantum-level surrogate model.
Following pretraining, we conducted fine-tuning on the SPICE dataset, which is characterized
by higher precision and presents more significant modeling challenges. Compared to nablaDFT,
SPICE offers significantly enhanced molecular diversity and quantum accuracy. We excluded ionic
species and non-neutral molecules from the SPICE dataset, ultimately retaining approximately 1.9
million conformations. These were split into a 95:5 train-test ratio, with the training set comprising
around 1.8 million conformations, covering molecules with 2 to 110 atoms. In addition to its
broader chemical diversity, SPICE includes dimers, solvated biomolecular conformations, and
peptides, making it more representative of real-world MD scenarios (as shown in Supplementary
Table 2). On the SPICE test set, LiTEN-FF achieved MAEs of 0.57 milliHatree for energy and 0.48
milliHatree/Å for force, corresponding to approximately 0.358 kcal/mol and 0.301 kcal/mol/Å,
respectively. Although these errors are slightly higher than those observed on the nablaDFT test set,
they still fall within the DFT-level error range, underscoring the model’s ability to maintain high
predictive accuracy even under distributional shifts. It is noteworthy that all evaluations in this study
were conducted using a molecule-level error metric, i.e., MAE computed over the entire molecular
energy and forces. This assessment strategy better reflects the model’s real-world performance in
applications such as drug-like molecule modeling, conformer optimization, and MD simulations.
During the SPICE finetuning phase, we particularly focused on evaluating the model’s ability to
adapt to variations in molecular scaffolds, conformational perturbations, and changes in electronic
environments. Through this fine-grained optimization on a high-accuracy dataset, we further
enhanced the model’s predictive accuracy and robustness in real drug discovery tasks. This approach
effectively transitions the model from broad-coverage general learning to high-precision specialized
modeling.
Exploration of Molecular Dynamics Simulations in Vacuum Systems with LiTEN-
FF
In molecular simulations, efficiently and accurately obtaining molecular conformations, geometric
parameters, and torsional energy profiles is a fundamental capability that underpins core tasks in
drug design and materials research45. Conformer optimization is essential for identifying stable
molecular structures and serves as the foundation for downstream applications such as docking and
property prediction. This is particularly critical in high-throughput drug screening scenarios, where
both speed and accuracy are imperative46. During MD simulations, the monitoring of fluctuations
in bond lengths and angles offers valuable insight into a molecule’s thermodynamic stability and
structural flexibility, making it a key indicator for evaluating the physical fidelity of force fields.
Furthermore, torsion angle scanning enables detailed a detailed characterization of intramolecular
energy variations, playing an indispensable role in conformational sampling, potential energy surface
construction, and force field parameterization. Therefore, focusing on these three critical tasks, we
systematically evaluated the applicability of the machine learning-based LiTEN-FF model under
vacuum conditions. The assessment encompasses two primary variants of the LiTEN-FF model,
namely LiTEN-FF Nab and LiTEN-FF SPICE, evaluated across three core tasks: conformer
optimization, geometric property analysis, and torsional energy profile prediction.
To systematically assess the applicability of LiTEN-FF Nab in the conformer optimization of
drug molecules, we conducted an in-depth performance analysis across multiple dimensions. A test
set of 300 molecules, with atom counts ranging from 8 to 62, was selected. These molecules were
excluded from the training set, thereby ensuring both the validity and representativeness of the
evaluation47. This size range covers the majority of commonly encountered drug-like compounds,
lending the study substantial practical relevance.
All initial conformers were uniformly generated using RDKit48, which eliminates biases arising
from different starting geometries and enables a fair comparison among various methods. The
optimization was performed via the ASE interface, with a uniform maximum force convergence
threshold of 0.05 eV/Å. This value represents the commonly used standard and can be adjusted as
needed. In terms of accuracy, the conformers optimized by LiTEN-FF Nab demonstrated excellent
agreement with the reference structures optimized by DFT (ωB97X-D/def2-SVP). The mean RMSD
between the two sets of structures was 0.05 Å, indicating the model’s ability to accurately reproduce
DFT-optimized geometries. More specifically, as shown in Fig. 3A, 31.2% of the molecules had
RMSD values below 0.01 Å, indicating near-exact reproduction of DFT geometries. Furthermore,
79.7% of the molecules had RMSD values below 0.05 Å, demonstrating very low structural deviation
for the majority of the test set. Additionally, 89.8% had RMSD values below 0.1 Å, further validating
the overall robustness and reliability of the optimized results.
The energy profiles of the optimization trajectories generated by LiTEN-FF Nab closely matched
those obtained from DFT, implying that the model not only accurately locates energy minima but
also captures the energy landscape along the optimization path, as shown in Supplementary Fig. 1.
This capability is essential for tasks such as conformer searching and MD simulations. Notably,
LiTEN-FF Nab demonstrated over a 1000-fold speedup compared to high-accuracy DFT methods.
This substantial boost in computational efficiency could greatly reduce the drug discovery timeline
and provides a strong technical foundation for high-throughput screening and virtual screening
pipelines. While the current evaluation is limited to molecules composed of elements included in
the training set, future work will aim to extend the model to systems containing metal atoms and
other complex functional groups. Incorporating solvent effects and temperature dependence into
the optimization process also represents a promising direction for exploring the model’s stability
under more realistic conditions. Furthermore, integrating LiTEN-FF Nab with MD simulations to
evaluate its performance in capturing dynamic conformational changes over extended timescales
could advance its applications in drug design and materials science.
To further assess LiTEN-FF’s ability to model dynamic geometric behavior, we conducted MD
simulations under vacuum conditions on 18 organic molecules not present in the training data.
These molecules represent a range of sizes (≤ 25 atoms) and contain typical organic elements such
as C, H, O, N, F, Cl, Br, and S. The simulations, lasting between 50 to 150 ps per molecule, were
designed to test the model’s fidelity in reproducing time-evolving molecular geometries and internal
coordinate distributions.
A comparative statistical analysis was conducted between the LiTEN-FF-predicted geometry
distributions and those derived from the DFT-based MD simulations. Specifically, we computed
Kullback-Leibler (KL) divergence metrics for bond lengths and bond angles across all molecules. As
shown in Fig. 3B-3C, the KL divergence remained below 0.01 in nearly all cases, indicating that
LiTEN-FF closely reproduces the geometric ensembles generated by DFT. Bond length distributions
showed tight peaks around equilibrium values, and angular distributions captured both symmetric
and asymmetric features of molecular flexibility. Furthermore, LiTEN-FF was able to maintain
overall molecular stability and structural integrity across simulation timescales, without introducing
artificial bond distortions or unphysical behaviors, an issue sometimes observed in less robust ML
potentials. We also present molecule-specific statistical summaries (Supplementary Table 3) and
visualizations (Supplementary Fig. 2-3) of the geometric property distributions, including variance,
skewness, and kurtosis. These analyses reveal consistent trends: smaller molecules typically exhibit
sharper and more narrowly clustered distributions, while slightly larger or more flexible molecules
show broader but still DFT-consistent behavior. This suggests that LiTEN-FF not only accurately
captures equilibrium geometries but also reflects the intrinsic thermal fluctuations in molecular
conformations, making it suitable for dynamic simulations that require physical realism over
extended timescales.
TorsionNet20632 is a high-accuracy molecular benchmark designed to evaluate the prediction
and ranking of dihedral energy profiles. The benchmark comprises 206 subsets, each containing 24
conformers, with energies computed at the CCSD(PT)/def2-TZVP level of theory. Performance
comparison on the the benchmark reveals key insights into the ability of neural network potential
models to accurately reproduce torsional energy profiles, a task that plays a vital role in drug design.
Traditional density functional theory methods, such as ωB97M D3BJ combined with the def2
TZVPPD basis set, remain the most accurate, achieving a MAE of 0.15 kcal per mol and serving as
the reference standard. Nevertheless, recent neural network models such as Egret One and the
MACE-OFF series have shown excellent performance, with MAEs ranging from 0.19 to 0.24 kcal
per mol and correlation coefficients exceeding 0.98 for both R squared and Spearman metrics. These
results greatly reduce the accuracy gap with DFT while delivering significant improvements in
computational efficiency.
Notably, as shown in Table 3, our LiTEN-FF SPICE model achieved an MAE of 0.19 kcal/mol
on this task, outperforming not only semi empirical approaches such as GFN2 xTB but also
traditional DFT methods like B3LYP-D3BJ/6-31G(d). This highlights LiTEN-FF SPICE’s ability to
effectively capture torsional energy landscapes shaped by complex electronic and steric interactions.
The pretrained LiTEN-FF Nab model, while slightly less accurate with an MAE of 0.54 kcal/mol,
still outperforms many earlier neural models. This highlights the critical role of training data quality,
as the accuracy of the model is closely tied to that of its training targets. The substantial performance
improvement observed when LiTEN-FF Nab is fine-tuned to LiTEN-FF SPICE clearly illustrates the
importance of domain-specific fine-tuning.
Importantly, the consistently high correlation metrics across different LiTEN-FF variants
demonstrate that the models not only provide accurate absolute energy predictions but also preserve
the correct ranking of torsional conformers. This is essential for downstream applications such as
conformer selection and force field parameterization.
In conclusion, LiTEN-FF SPICE delivers DFT-level accuracy in high-throughput torsion scan
tasks while drastically reducing computational costs. Compared to models like Egret-1 and MACE-
OFF, LiTEN-FF also offers a significant speed advantage since it does not rely on spherical harmonic
tensor computations. These results position LiTEN-FF SPICE as a practical and efficient alternative
for torsional energy modeling in real-world molecular design workflows.
LiTEN-FF MD Speed and Periodic Water Box Simulation
To comprehensively evaluate the computational efficiency of LiTEN-FF, we conducted a systematic
benchmark study against several widely used neural force field models, including MACE-OFF(L),
MACE-OFF(M), DPA2, and AIMNET2, across molecular system containing 100 to 1500 atoms. All
evaluated on the same NVIDIA GeForce RTX 4090 GPU. As shown in Fig. 4A, LiTEN-FF
consistently demonstrates the lowest computational overhead across all system sizes. For smaller
systems containing up to 300 atoms, it completes 10,000 MD steps in less than 200 seconds,
outperforming MACE-OFF23(L), which starts with a runtime of approximately 100 seconds but
rises sharply to over 500 seconds as the system size exceeds 300 atoms. Similarly, MACE-OFF23(M)
exceeds 400 seconds beyond this scale. For larger systems, such as those with 1200 atoms, LiTEN-
FF maintains a runtime below 400 seconds, achieving up to a ten-fold speedup compared to MACE-
OFF23(L), which exceeds 3500 seconds for the same system size.
This performance gap becomes increasingly evident as the system size expands, highlighting the
scalability advantages of LiTEN-FF. While other models tend to experience a rapid increase in
computational cost as molecular size grows, LiTEN-FF maintains efficient inference through its TQA
mechanism, which compactly encodes higher-order geometric correlations with linear
computational complexity. This architectural efficiency allows LiTEN-FF to scale gracefully, offering
substantially faster runtimes without sacrificing predictive accuracy.
Beyond speed, LiTEN-FF achieves efficiency without compromising physical fidelity. To validate
this, we evaluated LiTEN-FF’s performance in periodic systems using a cubic water box consisting of
348 water molecules. Given water’s critical role as both a fundamental molecular system and a
benchmark for evaluating biomolecular force fields, accurate modeling under periodic boundary
conditions is essential. To this end, we explored finetuning LiTEN-FF on the high-fidelity dataset
from Chen et al.49, which captures accurate thermodynamic behavior across diverse molecular
configurations. We then conducted 10 ns MD simulations to calculate the radial distribution
function (RDF) of oxygen-oxygen distances (Fig. 4B). The LiTEN-FF model finetuned on the ab
initio water dataset yielded RDFs that closely aligned with experimental measurements50, accurately
capturing both the first RDF peak and the subsequent hydration shells. In contrast, models trained
solely on the SPICE dataset, such as LiTEN-FF (SPICE) and MACE-OFF23(M), demonstrated nearly
identical deviations, including an overestimation of the first hydration valley (3-4 Å) and an
underestimation of the first RDF peak (around 4.5 Å). These systematic discrepancies are primarily
attributed to the limitations of the SPICE dataset, particularly its insufficient sampling of the water
configurational space and limited thermodynamic accuracy.
In addition to the RDF, we further evaluated whether LiTEN can accurately reproduce key
thermodynamic and dynamical properties of liquid water, including the density, heat of vaporization,
and self-diffusion coefficient. As shown in Supplementary Fig. 4, the densities obtained from LiTEN-
FF are in closer agreement with experimental measurements than those predicted by MACE-
OFF23(M) and MACE-OFF24(M) over the temperature range of 273–333 K. At 298 K, the heats of
vaporization obtained with LiTEN, MACE-OFF23(M), and MACE-OFF24(M) are 12.92 ± 0.04,
12.84 ± 0.14, and 12.88 ± 0.07 kcal·mol⁻¹, respectively (exp. 10.52 kcal·mol⁻¹). The MSD analysis
in Supplementary Fig. 5 compares the self-diffusion coefficients (𝐷) of water at 298 K computed by
LiTEN, MACE-OFF23(M), and MACE-OFF24(M) simulations, yielding D=0.240, 0.245, and 0.280
Å²·ps-1 (exp. 0.230 Å²·ps-1). This demonstrates that LiTEN not only reproduces the thermodynamic
properties of liquid water with high fidelity but also offers superior accuracy in describing water
diffusion, a property long recognized as particularly challenging to capture in molecular simulations.
Taken together, these results confirm that LiTEN-FF offers a compelling balance between
computational speed and physical fidelity, especially in large-scale and periodic systems. Its linear-
scaling architecture not only reduces simulation cost but also enables high-throughput exploration
of complex molecular environments, making it a strong candidate for next-generation machine
learning force fields in both vacuum and condensed-phase simulations.
Assessment of LiTEN-FF for Free Energy Surface Modeling in Vacuum and
Periodic Solvent Conditions
The free energy surface (FES) is a fundamental concept for characterizing the conformational
dynamics of biomolecules. It quantifies the relative stability of different conformations and maps
the kinetic pathways and energy barriers between them. Accurate reconstruction of the FES is
essential for elucidating mechanisms such as protein folding, enzymatic reactions, and ligand
binding. However, achieving both high accuracy and computational efficiency in FES calculations
remains a major challenge in molecular simulations and a critical benchmark for assessing force
fields and machine learning potentials. In this study, the performance of the LiTEN-FF model is
evaluated in the context of FES prediction, with a focus on its sampling efficiency and energy
estimation accuracy.
The alanine dipeptide (ala2) has long been regarded as a canonical system for testing
biomolecular simulation methods due to its well-characterized conformational states. Utilizing the
LiTEN-FF model, we performed 10 ns enhanced sampling simulations to characterize the backbone
dihedral angles 𝜑 and 𝜓, as shown in Fig. 5A and Supplementary Fig. 6. The implementation of
the on-the-fly probability enhanced sampling (OPES) technique accelerated the exploration of rarely
visited conformations, enabling the rapid convergence of free energy surface51. Across three
independent parallel trajectories, the free energy difference between key metastable states was
approximately 9.62 kJ/mol, with a statistical uncertainty significantly lower than 0.5 k T ,
b
underscoring the model’s excellent reproducibility and statistical robustness. Crucially, the free
energy profile predicted by LiTEN-FF aligns nearly perfectly with high-level DFT calculations,
therefore validating the model’s outstanding accuracy. This achievement demonstrates LiTEN-FF’s
ability to combine theoretical rigor with computational efficiency, providing a reliable foundation
for simulating more complex molecular systems.
Building on its proven performance in vacuum systems, LiTEN-FF was further applied to the
more challenging system of alanine tetrapeptide (ala4) in aqueous solution. This peptide exhibits
extensive conformational flexibility governed primarily by three backbone dihedral angles
𝜑 ,𝜑 ,𝑎𝑛𝑑 𝜑 51. We conducted 20 ns OPES simulations using these dihedrals as collective
1 2 3
variables, and our focus was on the slower-relaxing 𝜑 −𝜑 subspace to reveal the dominant
2 3
conformational basins. The free energy landscapes generated from LiTEN-FF (fine-tuned on the
SPICE dataset) and the MACE-OFF model show remarkable agreement, consistently identifying
four metastable states, as shown in Fig. 5B. The basin defined by 𝜑
2
< 0 and 𝜑
3
< 0 corresponds
to the lowest free energy minimum, and it features a structure akin to an antiparallel β-sheet, thereby
representing the peptide’s most stable conformation. Notably, LiTEN-FF predicts a slightly broader
basin compared to MACE-OFF, suggesting a more nuanced sampling of conformational
heterogeneity. Furthermore, LiTEN-FF exhibits superior computational efficiency relative to MACE-
OFF, highlighting its practical advantages for large-scale biomolecular simulations. These results
collectively confirm LiTEN-FF’s robustness and predictive power in accurately reproducing complex
biomolecular free energy landscapes within realistic solvent environments.
We performed unbiased NPT simulations using LiTEN-FF and MACE-OFF23(M) for both the
alanine tripeptide (Ala3, 20 ns) and alanine pentapeptide (Ala5, 10 ns), and calculated the ³J-
coupling constants from the Ramachandran distributions for direct comparison with experimental
data. The periodic water box system for Ala5 contained 6,869 atoms, while the Ala3 system
contained 1,026 atoms. For the Ala5 water box, simulations were carried out on a single NVIDIA
H100 GPU for 10 ns, taking 12 days for LiTEN-FF and 32 days for MACE-OFF23(M). Based on
our rough estimates, running MACE-OFF24(M) would require approximately 57 days, which was
computationally prohibitive; therefore, it was not included in the comparison. The results,
summarized in Table 4, demonstrate that both LiTEN-FF and MACE-OFF accurately reproduce the
3J-coupling constants. Moreover, as discussed above, LiTEN achieves comparable accuracy with
significantly higher computational efficiency, further highlighting the advantages of our approach.
In addition, Supplementary Fig. 7-10 present the Ramachandran distributions and the time
evolution of the RMSD of the trajectories with respect to the first frame for both LiTEN-FF and
MACE-OFF23. These results clearly show that LiTEN exhibits a noticeably faster conformational
diffusion compared to MACE-OFF23, further emphasizing its superior sampling efficiency.
Conformer search with LiTEN-FF
Generating multiple conformers of drug-like molecules is essential in drug design, as it enables
accurate modeling of ligand-target interactions by capturing the molecular flexibility critical for
binding and biological activity. After demonstrating the reliability of LiTEN-FF in molecular
simulations, we applied it to the task of conformer generation. Specifically, we performed multiple
rounds of high-temperature and annealing simulations on an entire batch of diverse molecules to
explore the conformational space. After each annealing cycle, we optimized the resulting structures
using the FIRE algorithm52 to obtain well-relaxed conformers. The initial coordinates for a new
simulation round were taken from the last frame of the previous high-temperature simulation.
To determine whether the generation process had converged, we introduced a convergence
criterion on RMSD analysis. For each newly generated conformer, we computed the RMSD against
all previously generated conformers of the same molecule. If, over a certain number of consecutive
rounds, every newly generated conformer could be matched to an existing conformer within a
predefined RMSD threshold, we considered the generation to have converged and terminated
further simulations. This ensured that no substantially new conformations were being sampled,
thereby indicating sufficient coverage of the conformational space. After the simulations were
completed, we removed redundant conformers based on the predefined RMSD threshold.
To demonstrate the efficacy of our method in thoroughly exploring the conformational space,
we performed t-SNE dimensionality reduction on the conformers generated for a benchmark set of
50 molecules and visualized the results using t-SNE plots (Fig. 6A). These conformations not only
exhibit structural diversity but also span across distinct energy basins, indicating broad coverage of
the underlying energy landscape. The widespread distribution of points in the t-SNE plot highlights
the method’s strength in capturing both high- and low-energy conformers, reflecting its ability to
overcome energy barriers and explore multiple metastable states. For single-molecule tasks, we
evaluated three representative drug molecules: Amoxicillin, Cloxacillin, and Linezolid, and
compared our results with those obtained using the CREST level of theory. As illustrated in Fig. 6B,
our method generates a greater number of diverse conformers than CREST, capturing a broader
range of structure variations. Additionally, the conformers produced by our approach exhibit more
detailed sampling, as reflected by lower RMSD values and better configurational convergence. In
terms of computational efficiency, our method is markedly faster, requiring less than half the
runtime of CREST across all three molecules (Fig. 6C).
We evaluated the efficiency of our conformer generation workflow using 1, 10, 30, 50, and 100
drug-like molecules with distinct structures. Each test involved 250 simulation rounds, with each
round consisting of 100 steps of high-temperature dynamics followed by 100 steps of annealing.
Additionally, the RMSD threshold for the convergence was detected in five consecutive rounds.
Notably, the inference speed of LiTEN-FF does not scale linearly with the number of atomic
coordinates, making it particularly well-suited for large-batch conformer generation. The detailed
workflow is illustrated in Fig. 6D. As shown in Fig. 6E, increasing the batch size to 100 accelerates
the conformer generation process by approximately 10-fold compared to generating conformers one-
by-one. All benchmarks were performed on an NVIDIA GeForce RTX 4090 GPU with 24 GB of
memory.
In summary, our generative method is capable of rapidly generating a large number of reliable
conformations for a wide range of drug-like molecules, significantly advancing the efficiency and
effectiveness of conformational sampling in drug discovery.

Methods
Invariance of Energy and Equivariance of Forces
In the fields of MD and geometric deep learning, it is crucial to ensure that the physical quantities
predicted by models satisfy fundamental symmetry requirements53. This section provides rigorous
mathematical proofs demonstrating the fundamental symmetry properties of the key physical
quantities in our framework.
The scalar energy function exhibits rotational invariance, satisfying:
E(𝑅𝒖) = E(𝒖). (8)
for any rotation matrix 𝑅. This invariance stems from two crucial properties:
Vector norm preservation: |𝑅𝒖 | = |𝒖 |, Angle cosine preservation:
𝒊 𝒊
𝑅𝒖 ⋅𝑅𝒖
cos𝜃 𝑗 ′ 𝑖𝑘 = 𝒋𝒊 𝒌𝒊 = cos𝜃 𝑗𝑖𝑘 . (9)
|𝑅𝒖 ||𝑅𝒖 |
𝒋𝒊 𝒌𝒊
The cross products term maintains rotational invariance:
(𝒖 ×𝒖 )⋅(𝒖 ×𝒖 ). (10)
𝒋 𝒋𝒌 𝒋𝒌 𝒌
Through the transformation property(𝑅𝑎)×(𝑅𝑏) = 𝑅(𝑎×𝑏) and orthogonality 𝑅ᵀ𝑅 = 𝐼,
we establish:
[(𝑅𝒖 )×(𝑅𝒖 )]·[( 𝑅𝒖 )×(𝑅𝒖 )] = (𝒖 ×𝒖 )·(𝒖 ×𝒖 ). (11)
𝒋 𝒋𝒌 𝒋𝒌 𝒌 𝒋 𝒋𝒌 𝒋𝒌 𝒌
The cubic norm expression inherits rotational invariance from its constituent invariant terms:
|𝒖 |3 = (∑cos𝜃 )·|𝒖 |. (12)
𝒊 𝑗𝑖𝑘 𝒊
Force equivariance follows from energy invariance through gradient transformation:
F(𝑅𝒖) = −𝛻 𝐸(𝑅𝒖) = −𝑅𝛻 𝐸(𝒖) = 𝑅F(𝒖), (13)
𝑅𝒖 𝑢
where we use the covariant gradient transformation 𝛻 = 𝑅𝛻 . These mathematical guarantees
𝑅𝑢 𝑢
ensure physical consistency under rotational transformations, which can be empirically verified by
applying random rotations to input coordinates and testing the invariance/equivariance conditions.
This theoretical foundation supports the development of reliable geometric learning models for
molecular simulations.
Tensorized Quadrangle Attention
In the Model Architecture section, we proposed an efficient and physically grounded mechanism
to capture dihedral interactions in molecular graphs, termed Tensorized Quadrangle Attention
(TQA). This method achieves linear complexity while maintaining high expressiveness by implicitly
modeling four-body geometric relationships during edge feature updates. The TQA mechanism is
formulated as:
𝐸𝑙+1 = [(𝒖⃗⃗⃗ ×𝒖⃗⃗⃗⃗⃗ )⋅(𝒖⃗⃗⃗⃗⃗ ×𝒖⃗⃗⃗ )]⊙SiLU(𝑊𝐸𝑙 +𝑏 ), (14)
𝑖𝑗 𝒊 𝒊𝒋 𝒊𝒋 𝒋 𝑒 𝑖𝑗 𝑒
where 𝑢⃗⃗⃗ ∈ 𝑅𝐶×𝟛 and 𝑢⃗⃗⃗ ∈ 𝑅𝐶×𝟛 are the directional vector embeddings of nodes 𝑖 and 𝑗, and
𝑖 𝑗
𝑢⃗⃗⃗⃗
𝑖
⃗
𝑗
∈ 𝑅𝐶×𝟛 is the normalized direction vector from node 𝑖 to node 𝑗 . The term(𝑢⃗⃗⃗
𝑖
×𝑢⃗⃗⃗⃗
𝑖
⃗
𝑗
)⋅
(𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗ ) captures a torsion-like relationship, analogous to a dihedral angle, but implemented in
𝑖𝑗 𝑗
a differentiable and vectorized form. This quantity modulates the updated edge features through a
Hadamard (element-wise) product with a nonlinear transformation of the previous edge state.
This formulation allows us to encode localized quadrangle interactions across edges, implicitly
modeling the torsional coupling between bonded atoms without the need for an explicit
enumeration of four-body terms. Since each message-passing step already includes information from
both node directions and edge vectors, the effective receptive field is enlarged. With a cutoff radius
of 5 Å for edge construction, this tensorized attention spans roughly 10 Å in a single layer, enabling
the network to capture mid- to long-range structural correlations efficiently.
To incorporate node-level information and further enrich the representation, we apply a learned
attention mechanism. Specifically, attention scores are computed as:
𝑎𝑙 = SiLU (ℎ𝑙 +𝐸𝑙 +ℎ𝑙)⊙ α, (15)
𝑖𝑗 𝑖 𝑖𝑗 𝑗
where ℎ𝑙 ∈ 𝑅𝐶 and ℎ𝑙 ∈ 𝑅𝐶 are the scalar features of atoms 𝑖 and 𝑗 , and α ∈ 𝑅𝐶 is a
𝑖 𝑗
trainable parameter vector that modulates attention weights. This attention 𝑎𝑙 ∈ 𝑅𝐶 is then gated
𝑖𝑗
by a continuous cutoff function based on the interatomic distance:
𝑎𝑙 = 𝑎𝑙 ⊙Radial_function(𝑑 ). (16)
𝑖𝑗 𝑖𝑗 𝑖𝑗
This hybrid design enables adaptive feature weighting that respects both learned chemical
relevance and geometric locality. The updated node and vector features are then computed as:
sca𝑙 = ℎ𝑙 ⊙𝐸𝑙 ⊙𝑎𝑙 , vec𝒍 = 𝒖𝒍 ⊙ℎ𝑙 +𝒖𝒍 ⊙ℎ𝑙. (17)
𝑚𝑖𝑗 𝑗 𝑖𝑗 𝑖𝑗 𝒎𝒊𝒋 𝒋 𝑗 𝒊𝒋 𝑗
The expression sca𝑙 ∈ 𝑅𝐶represents scalar message passing, while vec𝑙 ∈ 𝑅𝐶represents
𝑚𝑖𝑗 𝑚𝑖𝑗
vector message passing. Finally, both scalar and vector features are aggregated across neighbors using
scatter-based summation and residual connections:
ℎ𝑙 = scatter_sum (sca𝑙 )+ℎ𝑙, 𝒖𝒍 = scatter_sum( vec𝒍 )+𝒖𝒍. (18)
𝑚𝑖𝑗 𝒎𝒊𝒋
The expression ℎ𝑙 ∈ 𝑅𝐶 represents the scalar features of the node at the 𝑙 -th layer , respectively.
Similarly, 𝑢𝑙 ∈ 𝑅𝐶×𝟛 represents the vector features of the node at the 𝑙 -th layer. This construction
results in a fully vectorized, geometry-aware message passing scheme that models complex
conformational dependencies in a scalable and differentiable manner. By capturing implicit four-
body interactions through local vector operations, Tensorized Quadrangle Attention enhances the
model's capacity to represent torsional energy landscapes and anisotropic molecular features, all
while maintaining efficiency suitable for large-scale systems.
Scalar-Vector Fusion
This component serves to further integrate the vector and scalar node representations derived from
the Message Passing and TensorCrossFusion modules. The fusion mechanism is designed to leverage
both three-body and higher-order geometric interaction cues, whose theoretical underpinnings have
been rigorously justified in the Model section. For each node 𝑖, the vector feature 𝑢𝑙 ∈ 𝑅𝐶×𝟛and
𝑖
scalar feature be ℎ𝑙 ∈ 𝑅𝐶 are processed to capture complex geometric correlations and reinforce
𝑖
the bidirectional exchange of information between vectorial and scalar modalities. To this end, the
vector representation is first projected into two separate subspaces:
[ 𝒖 𝒍(𝟏) , 𝒖 𝒍(𝟐) ] = 𝐿𝑖𝑛𝑒𝑎𝑟 ( 𝒖𝒍 ). (19)
𝒊 𝒊 𝑣𝑒𝑐 𝒊
The expressions 𝑢 𝑙(1) ,𝑢 𝑙(2) ∈ 𝑅𝐶×𝟛 represent the split node vectors in the latent space at
𝑖 𝑖
layer 𝑙. These subspace projections are used to construct two geometric interaction terms: a vector
trilinear term, which captures directional alignment (akin to three-body angular interaction), and a
vector multi-body term, which encodes magnitude-based modulation (analogous to multi-body
geometric scaling). These terms are formally defined as:
3
𝒍(𝟏) 𝒍(𝟐)
VecTri = ∑ 𝒖 ⊙ 𝒖 , (20)
:,𝒋,∶ :,𝒋,∶
𝑗=1
3
2
VecMuti = (√∑( 𝒖 𝒍(𝟐) ) +𝜀) , 𝜀 = 10−8, (21)
:,𝒋,: 
𝑗=1
where VecTri, VecMuti ∈ 𝑅𝐶. These vector-derived signals are used to modulate scalar features.
The scalar representation is projected as:
[ℎ𝑙(1),ℎ𝑙(2),ℎ𝑙(3)] = Linear (ℎ𝑙). (22)
𝑠𝑐𝑎
Then, we compute the update residues:
Scalar update residue:
ℎ𝑙+1 = (VecTri+VecMuti)⊙ℎ𝑙(1) +ℎ𝑙(2). (23)
Vector update residue:
𝒖𝒍+𝟏 = 𝒖𝒍(𝟏) ⋅ℎ𝑙(3), (24)
where the scalar ℎ 𝑙(3) ∈ 𝑅𝐶is broadcast across the spatial dimension of 𝑢 𝑙(1) .
𝑖 𝑖
These update residues are added back to the original features in subsequent stages, enabling
tightly coupled scalar-vector representations that are sensitive to both directional geometry and
feature semantics. All symbols appearing in the manuscript are defined in Supplementary Table 4.
Dataset Construction
NablaDFT is a high-quality quantum chemical dataset specifically curated for pharmaceutical
modeling. After integrating multiple of its subsets, we obtained approximately 1.9 million molecules
and 16 million conformations. Each conformation is annotated with key quantum chemical
properties, including total energy, atomic forces, Hamiltonian matrices, and overlap matrices,
calculated at the ωB97X-D/def2-SVP level of theory. The dataset covers a wide range of
representative drug-like molecules with atom counts ranging from 8 to 62, and includes eight
common pharmaceutical elements: H, C, N, O, S, Cl, F, and Br.
SPICE is a recently released, high-fidelity quantum chemical dataset designed to benchmark the
generalization ability of AI-based force fields for small-molecule modeling. It includes over 110,000
small molecules and ~2 million conformations, covering 17 elements: H, Li, B, C, N, O, F, Na, Mg,
Si, P, S, Cl, K, Ca, Br, and I. Each conformation is annotated with high-accuracy quantum properties
(energy, atomic forces, and electronic descriptors) computed at the ωB97M-D3/def2-TZVPP level
of theory.
Model Training and Evaluation Setup
In the model evaluation experiments, we trained the models from scratch on the rMD17, MD22,
and Chignolin datasets. The dataset splits, number of layers, and hidden dimensionality have been
described in detail in the Results and Discussion section. For the rMD17 dataset, we used a cutoff
radius of 5 Å and an initial learning rate of 0.0001. For MD22 and Chignolin, the cutoff was set to
4 Å and the initial learning rate to 0.0002. The batch size was set to 2 for both rMD17 and MD22,
and 4 for Chignolin. All training procedures employed mean absolute error (MAE) as the loss
function, with a loss weight of 1 for energy and 99 for forces. Early stopping was used to select the
best-performing checkpoint on the validation set, which was then applied to test on all remaining
molecules to generate the reported results.
For the foundation model pretraining, we employed a 6-layer equivariant neural network with
a hidden dimension of 256. The batch size was set to 16 during training on the NablaDFT dataset
and 8 during fine-tuning on the SPICE dataset. The learning rate was scaled based on the number
of GPUs used in parallel. All training procedures were executed on eight NVIDIA GeForce RTX
4090 GPUs. Specifically, the model was trained for 100 epochs on NablaDFT and subsequently fine-
tuned for 50 epochs on SPICE. A 95:5 split ratio was used for partitioning the training and
validation sets during pretraining. Experimental settings for MD simulations are detailed in
Supplementary Note 1. Method Details. Specifically, molecular dynamics were conducted using ASE
with DFT-level or machine-learning potentials, employing Langevin dynamics for temperature
control and a Monte Carlo barostat for pressure regulation. Enhanced sampling was further achieved
via the On-the-fly Probability Enhanced Sampling (OPES) algorithm to accelerate the exploration of
configurational space across preselected collective variables.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Su and co-workers present a new machine learning potential (MLP) architecture in this manuscript. The key innovation is a\ntensorized quadrangle attention (TQA) mechanism designed to capture three- and four-body interactions. While the concept\nof explicitly incorporating such higher-order interactions to improve MLP accuracy is not new (as seen in ANI, MACE, and\nothers), nor is the idea of introducing approximations to improve efficiency (as in ViSNet), the TQA presents a novel way of\nreplacing tensor operations with linear operations. This seemingly modest technical innovation results in marginal\nimprovements in accuracy but substantial gains in training and inference speed, as demonstrated across various molecular\nmodelling tasks, including energy prediction, torsional scans, conformer search, and free energy surface reconstruction.\nGiven the importance of MLFF development, I find the work suitable for publication. Below are several comments that may\nhelp strengthen the manuscript:\n(1) While the incorporation of dihedral-aware directional messaging is an interesting idea, the overall architectural\ninnovation should be properly described. The main design components—edge-fusion attention, directional encoding via\nvector operations, and message-passing through Hadamard products—closely resemble the prior ViSNet model. The\nprimary extension here is the quadrangle attention, which incorporates four-body correlations via cross products of\ndirectional vectors. This is a reasonable and potentially meaningful addition, but its advantage over the ViSNet architecture\nis not convincingly demonstrated. Please provide additional comparisons and discussion.\n(2) The ability to simulate bulk water is essential. Based on these simulations, the authors should report how well the model\nreproduces key water properties such as density, heat of vaporization, and diffusion constant, in addition to the reported\nRDF.\n(3) The main practical advantage of the proposed architecture lies in its computational efficiency, which could enable routine\nnanosecond-scale MD simulations of condensed-phase systems. The authors are encouraged to apply LiTEN-FF to\nselected benchmark systems where MACE has been tested. For instance, simulations of (Ala)₅ or longer alanine peptides in\nwater can be used to compute NMR J couplings, offering a direct comparison with other models. Moreover, it would be\nvaluable to report actual MD run times for such systems, comparing LiTEN-FF with MACE and ViSNet.\n(4) MACE-OFF has been used for free energy perturbation calculations. Can LiTEN-FF support alchemical transformations\nas well? Please comment on this potential capability—or better, demonstrate it.\nMinor points\n(5) Some notational and presentational issues further limit clarity. For example, the terms V, E, and S in the model schematic\nare not clearly defined and do not correspond directly to any equations in the main text. Readers unfamiliar with ViSNet may\nfind this difficult to follow.\n(6) Line 449: a optimization threshold of 0.05 eV/A cannot be considered as a stringent criterion\n(7) Sometimes (line 545 and 547 for example) MACE is not corrected spelled.\n(8) In figure 1C, Bach should be Batch?", "reviewer_2": "In this article, the authors describe a message-passing neural network architecture for predicting energies and forces of\natomistic configurations. Their primary advancement is a specific formulation for the encoding of three- and four-body\ninteraction terms. Additionally, they train a foundation model using their proposed architecture for drug-like molecules and\ncreate a fine-tuned variant a smaller, higher-accuracy set of bio-organic molecules. They demonstrate state-of-the-art\nperformance on a number of downstream tasks. The paper is well-written and will be of interest to researchers working on\nbiomolecular applications, as well as researchers focused on neural network design for molecular modeling.\nI was unable to test the accompanying code due to errors described below. This should be remedied prior to publication of\nthe manuscript.\nAdditionally, I have the following suggestions for improving the manuscript:\nIn Figure 1, part C, it is unclear why some downstream task names are bolded. Additionally, there is a misspelling of 'batch'\nas 'bach'. Finally, the phrasing of the names of the 'Other Downstream Tasks' is unusual.\nIn the caption for Figure 1, an image of the 'final symbol' should be used if at all possible.\nIn line 181, it is stated that the descriptor is constructed from a weighted aggregation, yet equation (1) indicates a uniform\nweighting. Please clarify whether the terms of this sum have distinct weightings.\nIn Table 1 in the heading, the names of Allegro and ACE appear to be cut-off or misspelled.\nIn lines 286-289, the description appears contradictory. It is stated that all models use 6 layers, except NequIP and MACE,\nbut then it is stated that 6 layers are used for NequIP and then that 9 layers are used for VisNet. Please clarify.\nThe information in Figure 2 would be easier to understand as a table or bar chart.\nFor the Chignolin dataset, reporting speed and memory metrics for the models during inference/deployment rather than or in\naddition to the metrics during training would be helpful.\nIn the Torsion Angle Scanning section, please say a bit about how the TorsionNet206 dataset was generated, specifically\nwhat level of theory was used in its computation.\nFor the MD Speed test (starting in line 538), please clarify if all trials were completed on the same hardware. Additionally,\nplease clarify which MACE-OFF model was used for Figure 4B.\nThe word 'of' is repeated twice in line 662.\n(Remarks on code availability)\nI have encountered the following issues while attempting to set up an environment for LiTEN and run the examples:\n- Some packages in requirement.txt are not on PyPI and thus could not be found with `pip install -r requirement.txt` as\nsuggested in the installation instructions. This problem was alleviated by instead running `pip install -r requirement.txt --\nextra-index-url https://download.pytorch.org/whl/cu118 -f https://data.pyg.org/whl/torch-2.4.0+cu118.html `\n- After this, I experienced multiple version conflicts (e.g., e3nn=0.5.6 and e3nn=0.4.4 were requested but e3nn 0.5.6 requires\ne3nn=0.4.4).\n- I tried forcing an install with the requested versions, but then encountered issues with compatibility with my platform (e.g.,\nthe pyg-lib/torch-sparse wheel requires glibc >= 2.29, but my system only supports 2.28).\nPlease fix the installation process and test it on multiple systems if possible. You may consider replacing the requirement.txt\nfile with an environment.yml file, which allows you to specify additional conda channels where your packages can be found.\nYou can also create platform-specific conda-lock files if it is important to preserve the environment precisely.\nAdditionally, please remove hard-coded local paths (e.g., /home/suqun/…) from the codebase and replace them with relative\npaths or configurable options. At minimum, please direct the user to update these paths prior to running the example scripts.\nFinally, I recommend to make the LiTEN package installable (see https://packaging.python.org/en/latest/tutorials/packaging-\nprojects/), rather than instruct the user to add the directory to their path (e.g., in LiTEN/scripts/train_md17.pg, line 18).", "reviewer_3": "Summary: In this paper, a novel and innovative force field is designed based on the interactions among the triplets and\nquadraplets of atoms. This force field balances speed and accuracy, and is shown to be performant in a wide range of tasks,\nfrom optimization to sampling.\nSignificance: This paper is likely to be well perceived by the diverse readership of this journal, as it reports an interesting\nresult on the frontier of force field research.\nStrength:\n- The paper is well written and well illustrated.\n- The force field is accurate and efficient.\nLimitations:\n- The novelty is somewhat limited, and the specific designs of the functional forms need a bit more explanation and\nrationalization.\nQuestions:\n- Have the authors considered characterizing the expressive power of this functional form, and relating it to higher order\nspherical harmonics, or more easily, the representation theorem of the O(3) or SO(3) group? It appears that these functional\nforms can be easily derived from dot products.\nSuggestions:\n- I recommend that the authors explore the connection between the designed functional forms and the higher-order and\ncoupling terms in class-II classical force fields.\n- The equations 15-20 are written somewhat unclearly with confusing notations. The readability would be greatly improved if\nclearer mathematical expressions are used."}, "author_response": "Response to Reviewer #1\nWe sincerely thank the reviewer for their careful reading of our manuscript and for\nproviding constructive and insightful comments. The suggestions have helped us clarify\nthe novelty of our work, improve the presentation, and strengthen the experimental\nvalidation. Below we provide detailed responses to each point.\nQuestion 1. While the incorporation of dihedral-aware directional messaging is an\ninteresting idea, the overall architectural innovation should be properly described. The\nmain design components—edge-fusion attention, directional encoding via vector\noperations, and message-passing through Hadamard products—closely resemble the prior\nViSNet model. The primary extension here is the quadrangle attention, which\nincorporates four-body correlations via cross products of directional vectors. This is a\nreasonable and potentially meaningful addition, but its advantage over the ViSNet\narchitecture is not convincingly demonstrated. Please provide additional comparisons\nand discussion.\nResponse: Thank you for your comment. We have further elaborated on the innovative\naspects of our model in the revised manuscript. Our approach is centered on node\nvectors and systematically extends from three-body to four-body interactions, further\ngeneralizing to many-body interactions while incorporating coupling relationships\namong different many-body terms. As you noted, this concept is already discussed in\nthe manuscript, and VisNet adopts a similar strategy. However, VisNet’s design is\nlimited to four-body interactions, whereas our four-body formulation is theoretically\nmore complete and capable of capturing higher-order geometric correlations with\ngreater accuracy. It is worth emphasizing that we deliberately maintain simplicity in the\nmessage-passing modules of our framework, employing lightweight structures such as\nimplicit attention mechanisms to avoid unnecessary computational complexity. On this\nfoundation, LiTEN is able to preserve high predictive accuracy. Overall, our model\nprovides a more comprehensive theoretical framework, demonstrates superior\npredictive performance, and achieves significantly higher computational efficiency.\nFrom a methodological perspective, VisNet computes dihedral angles using projections\nof vectors in the edge-perpendicular space, as shown in the following equations:\n𝑡⃗⃗⃗⃗ =Rej (𝑣⃗⃗⃗ ) = ∑ 𝑁𝑖 Rej (𝑢⃗⃗⃗⃗⃗⃗ ), (1)\n𝑖𝑗 u⃗⃗⃗⃗i j 𝑖 𝑚=1 u⃗⃗⃗⃗i j 𝑖𝑚\n𝑡⃗⃗⃗ =Rej (𝑣⃗⃗⃗ ) = ∑ 𝑁𝑖 Rej (𝑢⃗⃗⃗⃗⃗ ), (2)\n𝑗𝑖 u⃗⃗⃗⃗j i 𝑖 𝑛=1 u⃗⃗⃗⃗j i 𝑖𝑛\n𝑁𝑖 𝑁𝑖 𝑁𝑖 𝑁𝑖 (3)\n⟨𝑡⃗⃗⃗⃗ ,𝑡⃗⃗⃗ ⟩ = ∑ ∑⟨Rej (𝑢⃗⃗⃗⃗⃗⃗⃗ ),Rej (𝑢⃗⃗⃗⃗⃗⃗ )⟩ = ∑ ∑𝑐𝑜𝑠𝜓 .\n𝑖𝑗 𝑗𝑖 𝑢⃗⃗⃗⃗𝑖⃗⃗𝑗 𝑖𝑚 𝑢⃗⃗⃗⃗𝑗⃗⃗ 𝑖 𝑖𝑛 𝑚𝑖𝑗𝑛\n𝑚=1𝑛=1 𝑚=1𝑛=1\nHowever, equation (3) is incomplete, as it neglects the magnitudes of the projected\nvectors. The differences in vector norms are non-negligible and significantly affect the\nresults. In contrast, our four-body calculations use cross products, explicitly accounting\nfor vector magnitudes:\ncos𝜙 = (𝑢⃗⃗⃗⃗ 𝑖 ⃗ 𝑘 ⃗ ×𝑢⃗⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 )⋅(𝑢⃗⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 ×𝑢⃗⃗⃗⃗ 𝑗 ⃗⃗ 𝑙) , (4)\n|𝑢⃗⃗⃗⃗ 𝑖 ⃗ 𝑘 ⃗ ×𝑢⃗⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 |⋅|𝑢⃗⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 ×𝑢⃗⃗⃗⃗ 𝑗 ⃗⃗ 𝑙|\n(5)\n∑ 𝑠𝑖𝑛𝜃 𝑠𝑖𝑛𝜃 𝑐𝑜𝑠𝜙 = ∑(𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗⃗⃗ )⋅(𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗ )\n𝑘𝑖𝑗 𝑙𝑗𝑖 𝑖𝑘 𝑖𝑗 𝑖𝑗 𝑗𝑙\n𝑘,𝑙∈𝒩(𝑖,𝑗) 𝑘,𝑙\n= (𝑢⃗⃗⃗ ×𝑢⃗⃗⃗⃗⃗ )⋅(𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗ ).\n𝑖 𝑖𝑗 𝑖𝑗 𝑗\nBy using the normal vectors of edges and nodes to compute dihedral angles, the\nresults have clear physical meaning. The cross product of edge and normal vectors\ncorresponds directly to the sine of the angle between the node vector and the edge,\nwhich complements the cosine values obtained from three-body interactions. This\nformula allows a complete description of four-body interactions. Furthermore, we\nextend this framework to many-body interactions:\n3 3 (6)\n|𝑢⃗⃗⃗ |3 = (𝐴+𝐵)3/2 ≈ 𝐴3/2+ 𝐴1/2𝐵+ 𝐴−1/2𝐵2.\n𝑖 2 8\nIn terms of accuracy, we conducted comprehensive evaluations on rMD17, MD22,\nand Chignolin datasets. Note that VisNet used a 9-layer network on rMD17 and MD22,\nwhile our model consistently uses 6 layers. Despite the smaller network, LiTEN\noutperforms VisNet in most metrics. For instance, on MD22 (7 molecules, 14 tasks),\nLiTEN surpasses VisNet in 9/14 tasks (Table R1). On rMD17 (10 molecules, 20 tasks),\nLiTEN matches or exceeds VisNet in 14/20 tasks (Table R2).\nTable R1: Direct comparisons of LiTEN and VisNet on MD22 (Energy (E, kcal/mol)\nand force (F, kcal/mol/Å)).\nMolecule Ac-Ala3-NHMe DHA Stachyose AT-AT\nProperty Energy Force Energy Force Energy Force Energy Force\natoms 42 56 87 60\nLiTEN 0.0597 0.0763 0.0876 0.0554 0.0845 0.0617 0.1263 0.0915\nVisNet 0.0636 0.0803 0.0741 0.0598 0.0915 0.0879 0.0708 0.0812\nMolecule AT-AT-CG-CG Buckyball_Catcher Dw_Nanotube\nProperty Energy Force Energy Force Energy Force\natoms 118 148 370\nLiTEN 0.1743 0.149 0.344 0.1143 0.6484 0.2139\nVisNet 0.196 0.148 0.537 0.201 0.601 0.292\nTable R2. Direct comparisons of LiTEN and VisNet on rMD17 (Energy (E, meV) and\nforce (F, meV/Å)).\nMolecule Aspirin Azobenzene Benzene Ethanol Malonaldehyde\nProperty Energy Force Energy Force Energy Force Energy Force Energy Force\nLiTEN 1.8 6.5 0.5 2.2 0.03 0.2 0.3 2.4 0.6 4.2\nVisNet 1.9 6.6 0.7 2.5 0.03 0.2 0.3 2.3 0.6 3.9\nMolecule Naphthalene Paracetamol Salicylic acid Toluene Uracil\nProperty Energy Force Energy Force Energy Force Energy Force Energy Force\nLiTEN 0.2 1.1 1.1 4.8 0.7 4.1 0.2 1.2 0.3 2.4\nVisNet 0.2 1.3 1.1 4.5 0.7 3.4 0.3 1.1 0.3 2.1\nIn terms of speed and overall performance, we followed the VisNet testing protocol,\nkeeping a 6-layer network with hidden size 128 for fairness. We additionally evaluated\na 6-layer LiTEN with hidden size 256. As shown in Table R3, LiTEN_128 achieves\n~30% faster training speed on Chignolin compared to VisNet, reduces GPU memory\nusage by ~60%, and improves accuracy by 17%. LiTEN_256 has slightly higher\nmemory usage, similar speed to VisNet_128, and higher accuracy.\nTable R3: Comprehensive performance comparison of LiTEN and VisNet on\nChignolin dataset\nModel ViSNet_128 LiTEN_128 LiTEN_256\nForce MAE (kcal/mol/Å) 0.4949 0.4094 0.3641\nGPU Memory (MiB) 9620 3912 13949\nTrain time consumption (ms) 31.68 22.75 31.21\nInference time consumption (ms) 11.62 8.09 10.53\nIn summary, these results demonstrate the superiority of our method and strategy\nin terms of completeness, accuracy, and efficiency. We have made the relevant\nmodifications in the ‘Model Architecture’ and ‘Training Efficiency Evaluation on\nthe Chignolin Dataset’ section as follows.\n“This module significantly enhances the model’s sensitivity to molecular conformational\ndistortions, rigid fragments, and deviations from coplanarity, allowing it to more accurately capture\ncomplex molecular structures such as conjugated systems, helical chains, and protein backbones. It\nis worth noting that VisNet constructs dihedral angles using the vector rejection method. Although\nthis approach can reflect local geometric relationships to some extent, it does not explicitly normalize\nthe magnitude of the rejection vectors when aggregating neighboring vectors, which may weaken the\ngeometric interpretability of the resulting expressions. In contrast, our formulation based on cross-\nproduct tensors provides a more direct correspondence to real geometric relationships and endows\nthe descriptors with clearer physical meaning. Furthermore, our model extends the four-body\ninteractions to higher-order many-body interactions and incorporates coupling relationships among\ndifferent interaction orders, thereby maintaining consistent physical representation in more complex\nmolecular systems.\nIn the LiTEN framework, all geometric and physical descriptors are constructed through tensor\noperations within atomic neighborhoods, enabling the three-body and four-body modules to achieve\nlinear computational complexity O(N), which is significantly more efficient than traditional\ncombinatorial approaches. At the same time, we intentionally maintain simplicity in the design of\nthe message-passing modules, employing lightweight structures such as implicit attention\nmechanisms to avoid unnecessary complexity. On this basis, LiTEN maintains excellent predictive\naccuracy and numerical stability. Overall, LiTEN combines high computational efficiency with\nstrong generalization ability, ensuring robust applicability across diverse chemical environments\nwithout relying on any predefined atom types.”\nFigure 2. Comparison of force prediction performance and computational efficiency of various models on the\nChignolin dataset. Metrics include: A: mean absolute error (MAE, kcal/mol/Å); B: GPU memory usage (MiB);\nC: training time (ms); D: inference time (ms). For MAE, lower values indicate higher prediction accuracy, while\nfor training/inference time and memory usage, lower values indicate greater computational efficiency.”\nQuestion 2. The ability to simulate bulk water is essential. Based on these simulations,\nthe authors should report how well the model reproduces key water properties such as\ndensity, heat of vaporization, and diffusion constant, in addition to the reported RDF.\nResponse: We thank the reviewer for this valuable suggestion. Beyond the radial\ndistribution function (RDF), other fundamental physical properties are equally crucial\nfor validating the accuracy of a force field. In response to the reviewer's comments,\nwe have further supplemented and systematically analyzed key thermodynamic and\ndynamic properties of liquid water—including density, heat of vaporization, and self-\ndiffusion coefficient—using MACE-OFF23(M) and MACE-OFF24(M) as reference\nmodels to comprehensively evaluate the performance of LiTEN-FF.\nAs shown in Supplementary Figure 4, across the temperature range of 273–333 K,\nthe water density predicted by LiTEN-FF exhibits significantly better agreement with\nexperimental measurements than MACE-OFF23 and MACE-OFF24, maintaining\nhigher consistency throughout the entire temperature range. At 298 K, the heat of\nvaporization calculated by LiTEN-FF, MACE-OFF23, and MACE-OFF24 are 12.92 ±\n0.04, 12.84 ± 0.14, and 12.88 ± 0.07 kcal·mol⁻¹, respectively (experimental value: 10.52\nkcal·mol⁻¹). Furthermore, the mean square displacement (MSD) analysis in\nSupplementary Figure 5 shows that at 298 K, the self-diffusion coefficients 𝐷 predicted\nby LiTEN-FF, MACE-OFF23, and MACE-OFF24 are 0.240, 0.245, and 0.280 Å²·ps⁻¹,\nrespectively (experimental value: 0.230 Å²·ps⁻¹). LiTEN-FF demonstrates closer\nagreement with experiment, successfully avoiding the over-diffusion commonly\nobserved in other models.\nOverall, LiTEN-FF excels in reproducing both thermodynamic and dynamic\nproperties of liquid water. It not only faithfully captures the experimental physical\ncharacteristics of liquid water but also demonstrates superior accuracy and stability in\nreproducing one of the most challenging properties in molecular simulations—the\ndiffusion of water molecules—highlighting its advantages in molecular force field\nmodeling.\nThe detailed description has been added to the section “LiTEN-FF MD Speed and\nPeriodic Water Box Simulation”, with the corresponding figures presented in\nSupplementary Figure 4 and Supplementary Figure 5.\n“In addition to the RDF, we further evaluated whether LiTEN can accurately reproduce key\nthermodynamic and dynamical properties of liquid water, including the density, heat of vaporization,\nand self-diffusion coefficient. As shown in Supplementary Figure 4, the densities obtained from\nLiTEN-FF are in closer agreement with experimental measurements than those predicted by MACE-\nOFF23 and MACE-OFF24 over the temperature range of 273–333 K. At 298 K, the heats of\nvaporization obtained with LiTEN, MACE-OFF23, and MACE-OFF24 are 12.92 ± 0.04,\n12.84 ± 0.14, and 12.88 ± 0.07 kcal·mol⁻¹, respectively (exp. 10.52 kcal·mol⁻¹). The MSD\nanalysis in Supplementary Figure 5 compares the self-diffusion coefficients (𝐷) of water at 298 K\ncomputed by LiTen, MACE-OFF23, and MACE-OFF24 simulations, yielding D=0.240, 0.245,\nand 0.280 Å²·ps-1 (exp. 0.230 Å²·ps-1). This demonstrates that LiTEN not only reproduces the\nthermodynamic properties of liquid water with high fidelity but also offers superior accuracy in\ndescribing water diffusion, a property long recognized as particularly challenging to capture in\nmolecular simulations.”\nSupplementary Figure 4. Temperature dependence of water density, comparing results from\nLiTEN-FF, MACE-OFF23, and MACE-OFF24 models with experimental data.\nSupplementary Figure 5. Comparison of mean-squared displacement (MSD) curves for water\nat 298 K computed using LiTEN and MACE-OFF23 machine-learned force fields. The linear\nslope of the MSD provides the self-diffusion coefficient 𝐷.\nQuestion 3. The main practical advantage of the proposed architecture lies in its\ncomputational efficiency, which could enable routine nanosecond-scale MD simulations\nof condensed-phase systems. The authors are encouraged to apply LiTEN-FF to selected\nbenchmark systems where MACE has been tested. For instance, simulations of (Ala)₅ or\nlonger alanine peptides in water can be used to compute NMR J couplings, offering a\ndirect comparison with other models. Moreover, it would be valuable to report actual\nMD run times for such systems, comparing LiTEN-FF with MACE and ViSNet.\nResponse: We thank the reviewer for this valuable suggestion and agree on the\nimportance of this set of experiments. Accordingly, following the methodology\ndescribed in the MACE-OFF original work, we performed unbiased NPT simulations\nfor the Ala3 system (1,026 atoms) over 20 ns, and, as per your suggestion, for the Ala5\nsystem (6,869 atoms) over 10 ns. MACE-OFF23(M) was used as the reference model.\nVisNet was not included because a foundation model directly applicable to this type of\nmolecular dynamics simulation has not yet been publicly released. Although AI2BMD\nis capable of performing related simulations, it divides proteins into multiple dipeptide\nunits and combines VisNet with classical polarizable force fields, resulting in a complex\nworkflow that differs substantially from both MACE-OFF and LiTEN-FF, making\ndirect comparison difficult. In terms of computational speed, we have already provided\nclear benchmarks in the first part of our study, and MD simulations simply repeat the\nsame framework computations, so LiTEN-FF maintains a significant speed advantage.\nMACE-OFF is sufficiently representative for comparison.\nFor the Ala5 water box system, 10 ns MD simulations were conducted on a single\nNVIDIA H100 GPU, where LiTEN-FF required 12 days and MACE-OFF23(M)\nrequired 32 days. Based on rough estimates, MACE-OFF24(M) would require\napproximately 57 days and was therefore not included due to computational limitations.\nWe also calculated the corresponding ³J-coupling constants for direct comparison with\nexperimental data. As summarized in Table 4, both LiTEN-FF and MACE-OFF\naccurately reproduce the experimentally measured ³J-coupling constants. Moreover,\nLiTEN-FF demonstrates significantly higher computational efficiency while\nmaintaining comparable predictive accuracy. Further analysis of Ramachandran\ndistributions and RMSD trajectories indicates that LiTEN-FF exhibits faster\nconformational diffusion compared to MACE-OFF23, further validating the\nadvantages of our method in practical applications.\nThese updates are reflected in the manuscript under the subsection “Assessment of\nLiTEN-FF for Free-Energy Surface Modeling in Vacuum and Periodic Solvent\nConditions”, with a newly added sub-subsection titled “Calculation and Validation\nof ³J-Coupling Constants in Alanine Tripeptide and Pentapeptide”.\n“Calculation and Validation of 3J-Coupling Constants in Alanine Tripeptide and\nPentapeptide. We performed unbiased NPT simulations using LiTEN-FF and MACE-\nOFF23(M) for both the alanine tripeptide (Ala3, 20 ns) and alanine pentapeptide (Ala5, 10 ns),\nand calculated the ³J-coupling constants from the Ramachandran distributions for direct comparison\nwith experimental data. The periodic water box system for Ala5 contained 6,869 atoms, while the\nAla3 system contained 1,026 atoms. For the Ala5 water box, simulations were carried out on a\nsingle NVIDIA H100 GPU for 10 ns, taking 12 days for LiTEN-FF and 32 days for MACE-\nOFF23(M). Based on our rough estimates, running MACE-OFF24(M) would require\napproximately 57 days, which was computationally prohibitive; therefore, it was not included in\nthe comparison. The results, summarized in Table 4, demonstrate that both LiTEN-FF and MACE-\nOFF accurately reproduce the 3J-coupling constants. Moreover, as discussed above, LiTEN achieves\ncomparable accuracy with significantly higher computational efficiency, further highlighting the\nadvantages of our approach. In addition, Supplementary Figures 7–10 present the Ramachandran\ndistributions and the time evolution of the RMSD of the trajectories with respect to the first frame\nfor both LiTEN-FF and MACE-OFF23. These results clearly show that LiTEN exhibits a\nnoticeably faster conformational diffusion compared to MACE-OFF23, further emphasizing its\nsuperior sampling efficiency.”\nTable 4. Summary of 3J-Coupling Constants Computed from Dihedral Distributions.\nsystem Residue Method 3J(HN,H\nα\n) 3J(HN,C\nβ\n) 3J(HN,C′)\nAla3 2 LiTEN 5.70(0.22) 2.42(0.09) 0.78(0.08)\nMACE-OFF23 5.46(0.17) 2.74(0.05) 0.67(0.01)\nEXP 5.68 2.39 1.13\nAla5 2 LiTEN 6.27(0.28) - -\nMACE-OFF23 5.95(0.25) - -\nEXP 5.59 - -\nMACE-OFF23 6.13(0.34) - -\nEXP 5.74 - -\nMACE-OFF23 6.13(0.21) - -\nEXP 5.98 - -\nSupplementary Figure 7. Ramachandran plots of alanine tripeptide obtained from 20 ns NPT\nsimulations performed using LiTEN-FF and MACE-OFF.\nSupplementary Figure 8. Time evolution of the RMSD of alanine tripeptide obtained from 20\nns NPT simulations performed using LiTEN-FF and MACE-OFF.\nSupplementary Figure 9. Ramachandran plots of alanine pentapeptide obtained from 10 ns\nNPT simulations performed using LiTEN-FF and MACE-OFF.\nSupplementary Figure 10. Time evolution of the RMSD of alanine pentapeptide obtained from\n10 ns NPT simulations performed using LiTEN-FF and MACE-OFF.\nQuestion 4. MACE-OFF has been used for free energy perturbation calculations. Can\nLiTEN-FF support alchemical transformations as well? Please comment on this potential\ncapability—or better, demonstrate it.\nResponse: We thank the reviewer for the valuable suggestion. As a machine learning-\nbased force field, LiTEN-FF provides continuous and differentiable predictions of\nenergies and forces, and is therefore, in principle, capable of supporting alchemical\ntransformations such as free energy perturbation (FEP) or thermodynamic integration\n(TI). Related frameworks already exist, for example, the integration of the MACE-OFF\nand ANI series into FEP calculations1. These approaches typically follow a QM/MM-\nlike strategy, constructing MLFF/MM models to partition the total energy as shown in\nEquation (7):\n𝐸 = 𝐸 +𝐸 +𝐸 . (7)\ntotal ML MM ML-MM\nThe cross-term 𝐸 is given by:\nML-MM\n𝑞 𝑞\n𝑖 𝑗\n𝐸𝐸 = ∑ ∑\n𝑀𝐿−𝑀𝑀 |𝑅𝑀𝑀 −𝑅𝑀𝐿|\n𝑖∈𝑀𝑀𝑗∈𝑀𝐿 𝑖 𝑗\n(8)\n𝐴 𝐵\n+∑ ∑ ( − ).\n𝑖∈𝑀𝑀 𝑗∈𝑀𝐿 12 6\n|𝑅𝑀𝑀−𝑅𝑀𝐿| |𝑅𝑀𝑀−𝑅𝑀𝐿|\n𝑖 𝑗 𝑖 𝑗\nIn such schemes, the core region is modeled with MLFF, the surrounding\nenvironment with MM, and their interaction is captured by the cross-term. This\napproach is flexible and allows the MLFF component to be replaced as long as the\nelement coverage is adequate. Given that our training dataset is comparable to or more\nextensive than that of MACE-OFF, and that LiTEN-FF has demonstrated robust\nperformance across diverse molecular dynamics tasks, we have strong reason to believe\nit can accurately capture energy changes from initial to target states. Although we have\nnot yet performed alchemical free energy calculations with LiTEN-FF and thus cannot\nprovide direct examples, its high computational efficiency and stable energy predictions\nmake it particularly well-suited for large-scale or complex alchemical studies. We plan\nto explore this application in future work to fully assess LiTEN-FF’s potential in drug\ndesign and molecular optimization.\nMinor points\nQuestion 5. Some notational and presentational issues further limit clarity. For example,\nthe terms V, E, and S in the model schematic are not clearly defined and do not\ncorrespond directly to any equations in the main text. Readers unfamiliar with ViSNet\nmay find this difficult to follow.\nResponse: We thank the reviewer for pointing this out. The symbols 𝑉, 𝐸 and 𝑆 are\nindeed shown in the model schematic, but they were not explicitly defined in the main\ntext. In the revised manuscript, we have added clear definitions for these symbols and\nclarified their corresponding mathematical expressions both in the figure caption and\nin the main text, making it easier for readers—especially those unfamiliar with\nViSNet—to understand the model architecture.\nWe have revised the legend of Figure 1 in the manuscript and added a unified\nexplanation of all symbols used in the text in Supplementary Table 4.\n“(A) The core architecture of LiTEN, illustrating the geometric computations and their corresponding physical\ninterpretations in terms of multi-body interactions. Here, ‘S’ denotes node scalar features, ‘V’ denotes node vector\nfeatures, and ‘E’ denotes edge features. The operations are defined as follows: ⨀ represents the Hadamard product,\n⨂ the cross product, ⨁ summation, and ⦿ denotes the dot product.”\nSupplementary Table 4. Glossary of notations.\nNotation Description\n𝑖; 𝑗; 𝑘; 𝑙 Node (Atom) index in molecular graph (e.g., i, j for adjacent atoms)\n𝑁(𝑖) Neighborhood of node i (set of adjacent atoms around atom i)\n𝑙; 𝐿 Index of network layers; Total number of network layers\n𝑆 Node scalar embedding (encodes non-directional local info)\n𝑉 Node vector embedding (encodes directional local info)\n𝐸 Edge embedding between node i and j\n𝑖𝑗\nℎ Node scalar embedding\n𝑟⃗⃗⃗ Relative position vector between node i and j\n𝑖𝑗\n𝑢⃗⃗⃗⃗⃗ Unit direction vector from node i to j\n𝑖𝑗\n𝑢⃗⃗⃗ First-order equivariant vector of node i (encodes directional anisotropy)\n𝑖\n𝑑 Distance between node i and j\n𝑖𝑗\n⊙ Hadamard product (element-wise product of matrices/vectors)\n⨂ Cross product of two vectors\n⦿ Dot (inner) product of two vectors\n⨁ Summation of two vectors or scalars\n𝑆𝑐𝑎𝑡𝑡𝑒𝑟 𝑠𝑢𝑚 Scatter-sum reduce function (aggregates)\n3𝐽 Spin-spin coupling constant (NMR property index for dynamic structure validation)\n𝐴 Self-interaction term in 3rd-order tensor features\n𝐵 Cross-interaction term in 3rd-order tensor features\n𝑉𝑒𝑐𝑇𝑟𝑖 Vector trilinear term for scalar-vector fusion\n𝑉𝑒𝑐𝑀𝑢𝑙𝑡𝑖 Vector multi-body term for scalar-vector fusion\n𝑅𝑎𝑑𝑖𝑎𝑙 𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛 Distance cutoff function (modulates attention weights by 𝑑 )\n𝑖𝑗\n𝑅 Rotation matrix\n𝐼 Identity matrix (used in proving rotation matrix orthogonality)\nQuestion 6. Line 449: a optimization threshold of 0.05 eV/A cannot be considered as a\nstringent criterion\nResponse: Regarding the optimization threshold of 0.05 eV/Å mentioned in line 449, we\nacknowledge that this value cannot be considered a stringent criterion. In practical applications,\nusers will be allowed to set the convergence threshold according to their specific requirements,\nthereby balancing computational accuracy and efficiency. We have added a note in the revised\nmanuscript to indicate that this threshold can be adjusted for different use cases.\nWe have made the corresponding revisions in the “Conformation Optimization Using\nLiTEN-FF Nab” section.\n‘The optimization was performed via the ASE interface, with a uniform maximum force\nconvergence threshold of 0.05 eV/Å. This value represents the commonly used standard and can\nbe adjusted as needed.’\nQuestion 7. Sometimes (line 545 and 547 for example) MACE is not corrected spelled.\nResponse: Regarding the spelling of MACE (e.g., lines 545 and 547), we have carefully\nchecked the entire manuscript and have consistently corrected it to the proper spelling,\n“MACE.”\nQuestion 8. In figure 1C, Bach should be Batch?\nResponse: Regarding “Bach” in Figure 1C, thank you for pointing this out. We have\ncorrected it to “Batch.”\nResponse to Reviewer #2\nWe sincerely thank the reviewer for recognizing the potential impact of our work on\nbiomolecular applications and neural network design for molecular modeling. We also\nappreciate the detailed feedback regarding the presentation of downstream tasks, figure\nclarity, methodological descriptions, and code reproducibility. These constructive\ncomments have been invaluable for improving both the scientific content and the\nusability of our code. Below, we provide detailed point-by-point responses to each\ncomment.\nQuestion 1. In Figure 1, part C, it is unclear why some downstream task names are\nbolded. Additionally, there is a misspelling of 'batch' as 'bach'. Finally, the phrasing of\nthe names of the 'Other Downstream Tasks' is unusual.\nResponse: We thank the reviewer for the careful comments. In Figure 1C, the bolded\ndownstream task names were originally intended to emphasize representative task\ncategories; however, we realize that this was not clearly indicated in the text. In the\nrevised manuscript, we have applied bold formatting consistently to all text in panel C\nto maintain visual uniformity. Additionally, the typographical error “bach” has been\ncorrected to “batch.” We have also rephrased the task names under “Other Downstream\nTasks” to make their descriptions more natural and clearer.\nQuestion 2. In the caption for Figure 1, an image of the 'final symbol' should be used if\nat all possible.\nResponse: We thank the reviewer for the valuable suggestion. We understand that the reviewer\nwould like the “final symbol” in Figure 1’s caption to be visually represented, allowing readers\nto more intuitively relate it to the elements in the figure. In the revised manuscript, we have\nupdated the legend of Figure 1 to include a visual depiction of this symbol, making the\ncorrespondence between the figure and its description clearer and more intuitive.\nThe specific modifications can be found in the legend of Figure 1.\n“(A) The core architecture of LiTEN, illustrating the geometric computations and their corresponding physical\ninterpretations in terms of multi-body interactions. Here, ‘S’ denotes node scalar features, ‘V’ denotes node vector\nfeatures, and ‘E’ denotes edge features. The operations are defined as follows: ⨀ represents the Hadamard product,\n⨂ the cross product, ⨁ summation, and ⦿ denotes the dot product.”\nQuestion 3. In line 181, it is stated that the descriptor is constructed from a weighted\naggregation, yet equation (1) indicates a uniform weighting. Please clarify whether the\nterms of this sum have distinct weightings.\nResponse: We sincerely thank the reviewer for the careful and insightful comments.\nThe issues you raised are entirely valid. In the main text, we presented the formula\nprimarily to illustrate the geometric construction principles of the graph network within\nour model. In the actual implementation, we introduce scalar weights for each (𝑢 )\n𝑖𝑗\nduring aggregation, enabling a non-uniformly weighted aggregation process. The\nrelevant implementation details are provided in the Methods section of the\nSupplementary Information, where we also discuss the mapping between our\nframework and spherical harmonic convolutions:\n“In LiTEN, geometric descriptors—including first-order node vectors, second- and third-order scalar\ninvariants, and four-body dihedral descriptors—are explicitly constructed to satisfy rotational\nequivariance or invariance. Each descriptor can be rigorously mapped to an irreducible\nrepresentation of SO(3) and interpreted in the context of spherical harmonic convolutions, ensuring\nboth theoretical consistency and practical expressive power.\nThe first-order node vector associated with atom (i) is defined as:\n𝑟⃗⃗⃗ (1)\n𝑖𝑗\n𝑢⃗⃗⃗ = ∑ 𝑤 ⋅ ,\n𝑖 𝑖𝑗\n|𝑟⃗⃗⃗ |\n𝑖𝑗\n𝑗∈𝒩(𝑖)\nwhere 𝑟⃗⃗⃗ = 𝑟⃗⃗ −𝑟⃗⃗ is the displacement vector to a neighboring atom (j). Under a rotation 𝑅 ∈\n𝑖𝑗 𝑗 𝑖\n𝑆𝑂(3) , this vector transforms equivariantly as 𝑢⃗⃗⃗ → 𝑅𝑢⃗⃗⃗ , corresponding to the 𝑙 = 1\n𝑖 𝑖\nirreducible representation. This vector encodes directional information about the local environment\nof atom (i) and directly corresponds to the 𝑙 = 1 component in a spherical harmonic convolution\nof the form:\n𝑟⃗⃗⃗ (2)\n∑ 𝑌𝑚(𝑟 ̂ )𝑊𝑚𝑓 , 𝑟 ̂ = 𝑖𝑗 ,\n1 𝑖𝑗 1 𝑗 𝑖𝑗\n|𝑟⃗⃗⃗ |\n𝑖𝑗\n𝑗∈𝒩(𝑖)\nwhere 𝑌𝑚(𝜃,𝜙) are the spherical harmonics and 𝑊𝑚 are learnable weights. By using unit\n1 1\n̂\nvectors 𝑟 , LiTEN achieves a computationally efficient approximation of the spherical harmonic\n𝑖𝑗\nbasis without explicit evaluation of the 𝑌𝑚 functions.’\noff or misspelled.\nResponse: We thank the reviewer for the reminder. Due to typesetting and rendering issues,\nthe names “Allegro” and “ACE” in Table 1 were partially truncated. In the revised version, we\nhave corrected the formatting to ensure that the model names are displayed fully and clearly.\nThe specific modifications can be found in Table 1.\nTable 1. Mean absolute errors on the rMD17 dataset. Energy (E, meV) and force (F, meV/Å) errors for different\nmodels trained on 950 configurations, validated on 50, and evaluated on the held-out test set.\nLiTEN MACE Allegro VisNet BOTNet NequIP GemNet ACE ANI PaiNN\nEnergy 1.8 2.2 2.3 1.9 2.3 2.3 (T/-Q ) 6.1 16.6 6.9\nAspirin\nForce 6.5 6.6 7.3 6.6 8.5 8.2 9.5 17.9 40.6 16.1\nEnergy 0.5 1.2 1.2 0.7 0.7 0.7 - 3.6 15.9 -\nAzobenzene\nForce 2.2 3.0 2.6 2.5 3.3 2.9 - 10.9 35.4 -\nEnergy 0.03 0.4 0.3 0.03 0.03 0.04 - 0.04 3.3 -\nBenzene\nForce 0.2 0.3 0.2 0.2 0.3 0.3 0.5 0.5 10.0 -\nEnergy 0.3 0.4 0.4 0.3 0.4 0.4 - 1.2 2.5 2.7\nEthanol\nForce 2.4 2.1 2.1 2.3 3.2 2.8 3.6 7.3 13.4 10.0\nEnergy 0.6 0.8 0.6 0.6 0.8 0.8 - 1.7 4.6 3.9\nMalonaldehyde\nForce 4.2 4.1 3.6 3.9 5.8 5.1 6.6 11.1 24.5 13.8\nEnergy 0.2 0.5 0.2 0.2 0.2 0.9 - 0.9 11.3 5.1\nNaphthalene\nForce 1.1 1.6 0.9 1.3 1.8 1.3 1.9 5.1 29.2 3.6\nEnergy 1.1 1.3 1.5 1.1 1.3 1.4 - 4.0 11.5 -\nParacetamol\nForce 4.8 4.8 4.9 4.5 5.8 5.9 - 12.7 30.4 -\nEnergy 0.7 0.9 0.9 0.7 0.8 0.7 - 1.8 9.2 4.9\nSalicylic acid Force 4.1 3.1 2.9 3.4 4.3 4.0 5.3 9.3 29.7 9.1\nEnergy 0.2 0.5 0.4 0.3 0.3 0.3 - 1.1 7.7 4.2\nToluene\nForce 1.2 1.5 1.8 1.1 1.9 1.6 2.2 6.5 24.3 4.4\nEnergy 0.3 0.5 0.6 0.3 0.4 0.4 - 1.1 5.1 4.5\nUracil\nForce 2.4 2.1 1.8 2.1 3.2 3.1 3.8 6.6 21.4 6.1\nQuestion 5. In lines 286-289, the description appears contradictory. It is stated that all\nmodels use 6 layers, except NequIP and MACE, but then it is stated that 6 layers are used\nfor NequIP and then that 9 layers are used for VisNet. Please clarify.\nResponse: We thank the reviewer for the reminder. In the original text, our initial intention\nwas to highlight that MACE and Nequip are somewhat special, as they are based on spherical\nharmonic convolutions and thus have particular settings for layers and hidden sizes. However,\nwe agree that this may cause ambiguity. In fact, we have fully followed the results and\nexperimental settings reported in their original papers. Therefore, we have revised the text in\n‘Tests of model accuracy and generalization for quantitative properties’ section as follows.\n“In the experimental setup, we uniformly used the optimal hyperparameter configurations\nrecommended in the literature for each model to ensure fairness in comparison. Except for VisNet\nand MACE, which used their recommended numbers of layers (9 layers and 2 layers respectively),\nother models employed a 6-layer network architecture with a hidden dimension size of 256.”\nQuestion 6. The information in Figure 2 would be easier to understand as a table or bar\nchart.\nResponse: Thank you for your suggestion. We have revised Figure 2, presenting it as a\nbar chart with the exact values labeled, and have included the corresponding raw data\nin the supplementary files.\nThe updated Figure 2 in the manuscript is as follows:\nFigure 2. Comparison of force prediction performance and computational efficiency of various models on the\nChignolin dataset. Metrics include: A: mean absolute error (MAE, kcal/mol/Å); B: GPU memory usage (MiB);\nC: training time (ms); D: inference time (ms). For MAE, lower values indicate higher prediction accuracy, while\nfor training/inference time and memory usage, lower values indicate greater computational efficiency.\nQuestion 7. For the Chignolin dataset, reporting speed and memory metrics for the\nmodels during inference/deployment rather than or in addition to the metrics during\ntraining would be helpful.\nResponse: Thank you for your suggestion. We agree that it is very meaningful and have\nincluded a comparison of inference speed in Figure 2. Regarding GPU memory usage\nduring inference, we found it to be similar to that during training, because these models\ncompute forces via gradients. Therefore, whether in training or inference, gradient\ninformation and related parameters still need to be stored.\nThe specific modifications in the manuscript are shown in the ‘Response to Question\n6’.\nQuestion 8. In the Torsion Angle Scanning section, please say a bit about how the\nTorsionNet206 dataset was generated, specifically what level of theory was used in its\ncomputation.\nResponse: We thank the reviewer for the suggestion. TorsionNet206, developed by Xiao\nJ et al.2, is a high-accuracy molecular benchmark for evaluating the prediction and\nranking of dihedral energy profiles. The benchmark consists of 206 subsets, each\ncontaining 24 conformers, with energies computed at the CCSD(T)/def2-TZVP level\nof theory. We have incorporated this description into the manuscript in the section\n“Torsion Angle Scanning: TorsionNet206 Dataset.”\n“TorsionNet206 is a high-accuracy molecular benchmark designed to evaluate the prediction and\nranking of dihedral energy profiles. The benchmark comprises 206 subsets, each containing 24\nconformers, with energies computed at the CCSD(T)/def2-TZVP level of theory.”\nQuestion 9. For the MD Speed test (starting in line 538), please clarify if all trials were\ncompleted on the same hardware. Additionally, please clarify which MACE-OFF model\nwas used for Figure 4B.\nResponse: Thank you for your suggestion. We have clarified in the manuscript that all\ntests were conducted on an NVIDIA RTX 4090 GPU, and we have specified the version\nof MACE-OFF used (MACE-OFF23(M)). The corresponding revision can be found in\nthe “LiTEN-FF MD Speed and Periodic Water Box Simulation” section:\n“To comprehensively evaluate the computational efficiency of LiTEN-FF, we conducted a\nsystematic benchmark study against several widely used neural force field models, including MACE-\nOFF(L), MACE-OFF(M), DPA2, and AIMNET2, across molecular system containing 100 to\n1500 atoms. All evaluated on the same NVIDIA RTX 4090 GPU.”\nFigure 4. (A) Runtime comparison of LiTEN-FF and baseline models for 10,000 MD steps across systems with\nvarying atom counts. (B) Performance of LiTEN-FF in a periodic water box, as shown by the radial distribution\nfunction between oxygen atoms (O–O).\nQuestion 10. The word 'of' is repeated twice in line 662.\nResponse: Thank you for your careful observation. We have revised the manuscript\naccordingly and removed the redundant “of.”\nRemarks on code availability:\nI have encountered the following issues while attempting to set up an environment for\nLiTEN and run the examples:\n- Some packages in requirement.txt are not on PyPI and thus could not be found with\n`pip install -r requirement.txt` as suggested in the installation instructions. This problem\nwas alleviated by instead running `pip install -r requirement.txt -extra-index-url\nhttps://download.pytorch.org/whl/cu118 -f https://data.pyg.org/whl/torch-\n2.4.0+cu118.html `\n- After this, I experienced multiple version conflicts (e.g., e3nn=0.5.6 and e3nn=0.4.4\nwere requested but e3nn 0.5.6 requires e3nn=0.4.4).\n- I tried forcing an install with the requested versions, but then encountered issues with\ncompatibility with my platform (e.g., the pyg-lib/torch-sparse wheel requires glibc >= 2.29,\nbut my system only supports 2.28).\nPlease fix the installation process and test it on multiple systems if possible. You may\nconsider replacing the requirement.txt file with an environment.yml file, which allows\nyou to specify additional conda channels where your packages can be found. You can\nalso create platform-specific conda-lock files if it is important to preserve the environment\nprecisely.\nAdditionally, please remove hard-coded local paths (e.g., /home/suqun/…) from the\ncodebase and replace them with relative paths or configurable options. At minimum,\nplease direct the user to update these paths prior to running the example scripts.\nFinally, I recommend to make the LiTEN package installable (see\nhttps://packaging.python.org/en/latest/tutorials/packaging-projects/), rather than\ninstruct the user to add the directory to their path (e.g., in LiTEN/scripts/train_md17.pg,\nline 18).\nResponse: We sincerely thank the reviewer for the valuable feedback on the usability\nof the LiTEN code, and we apologize for the issues caused by the previous environment\nsetup. The earlier environment included numerous packages required for molecular\ndynamics simulations and other model runs, which may have inadvertently introduced\nversion conflicts and installation difficulties. In response to the installation and\nenvironment issues raised, we have comprehensively optimized both the code and\ndocumentation. We provide three options for environment setup: First, the most stable\nand recommended method. A complete environment.yml file and the corresponding\nConda package are available on Zenodo (https://zenodo.org/records/17149388). Users\ncan simply download the package to their envs directory, extract it, and activate the\nenvironment for immediate use. Second, we provide a convenient bash script. Users\ncan install all dependencies by simply running:‘bash liten_env.sh’. Third, we also offer\na complete environment.yml file, allowing users to directly create the required\nenvironment via ‘conda env create -f environment.yml’. The approach helps avoid\ncommon issues with missing PyPI packages or version conflicts that may occur when\nusing a requirements.txt file. Additionally, all hard-coded local paths have been\nreplaced with relative paths or configurable parameters, so users do not need to modify\nsource code when running example scripts. By simply downloading the data files from\n(https://zenodo.org/records/17142012) into the corresponding data directory, the scripts\ncan be executed directly. These improvements greatly simplify environment setup and\nenhance code portability and reproducibility, with detailed operational instructions\nprovided in the GitHub repository (https://github.com/lingcon01/LiTEN/tree/master)\nREADME.\nTo ensure environment stability and cross-platform compatibility, we have tested\nthe Conda environment on a variety of hardware, including NVIDIA 4090, H100, and\nH800, verifying version compatibility of core dependencies such as PyTorch, PyG, and\ne3nn, and resolving potential installation conflicts. Although LiTEN is not yet packaged\nas an installable Python package, users can still conveniently run all example scripts\nand training workflows using the provided environment.yml and relative path\nconfigurations, fully reproducing the experiments reported in the manuscript, and\nenabling further development and extension on this foundation.\nResponse to Reviewer #3\nWe sincerely thank the reviewer for their thoughtful evaluation and for recognizing the\naccuracy, efficiency, and broad applicability of our proposed force field. We also\nappreciate the insightful suggestions regarding the explanation of the functional forms,\ntheir connection to higher-order terms, and the mathematical clarity of equations. These\ncomments have been extremely helpful in improving both the rigor and the presentation\nof our manuscript. Below, we provide detailed point-by-point responses to each comment.\nQuestion 1. Have the authors considered characterizing the expressive power of this\nfunctional form, and relating it to higher order spherical harmonics, or more easily, the\nrepresentation theorem of the O(3) or SO(3) group? It appears that these functional\nforms can be easily derived from dot products.\nResponse: We thank the reviewer for the insightful comment regarding the expressive\npower of the functional forms used in LiTEN and their relation to higher-order spherical\nharmonics as well as the representation theory of the O(3)/SO(3) group. In LiTEN,\ngeometric descriptors—including first-order node vectors, second- and third-order\nscalar invariants, and four-body dihedral descriptors—are explicitly constructed to\nsatisfy rotational equivariance or invariance. Each descriptor can be rigorously mapped\nto an irreducible representation of SO(3) and interpreted in the context of spherical\nharmonic convolutions, ensuring both theoretical consistency and practical expressive\npower.\nThe first-order node vector associated with atom (i) is defined as:\n𝑟⃗⃗⃗ (7)\n𝑖𝑗\n𝑢⃗⃗⃗ = ∑ 𝑤 ⋅ ,\n𝑖 𝑖𝑗\n|𝑟⃗⃗⃗ |\n𝑖𝑗\n𝑗∈𝒩(𝑖)\nwhere 𝑟⃗⃗⃗ = 𝑟⃗⃗ −𝑟⃗⃗ is the displacement vector to a neighboring atom (j). Under a\n𝑖𝑗 𝑗 𝑖\nrotation 𝑅 ∈ 𝑆𝑂(3), this vector transforms equivariantly as 𝑢⃗⃗⃗ → 𝑅𝑢⃗⃗⃗ , corresponding\n𝑖 𝑖\nto the 𝑙 = 1 irreducible representation. This vector encodes directional information\nabout the local environment of atom (i) and directly corresponds to the 𝑙 =\n1 component in a spherical harmonic convolution of the form:\n∑ 𝑌𝑚(𝑟 ̂ )𝑊𝑚𝑓 , 𝑟 ̂ = 𝑟⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 , (8)\n𝑗∈𝒩(𝑖) 1 𝑖𝑗 1 𝑗 𝑖𝑗\n|𝑟⃗⃗⃗ 𝑖 ⃗⃗ 𝑗 |\nwhere 𝑌𝑚(𝜃,𝜙) are the spherical harmonics and 𝑊𝑚 are learnable weights. By\n1 1\n̂\nusing unit vectors 𝑟 , LiTEN achieves a computationally efficient approximation of\n𝑖𝑗\nthe spherical harmonic basis without explicit evaluation of the 𝑌𝑚 functions3, 4.\nvectors:\n|𝑢⃗⃗⃗ |2 = ∑ 𝑢⃗⃗⃗⃗⃗ ⋅𝑢⃗⃗⃗⃗⃗ = ∑ cos𝜃 , (9)\n𝑖 𝑗,𝑘∈𝒩(𝑖) 𝑖𝑗 𝑖𝑘 𝑗,𝑘 𝑗𝑖𝑘\nwhere 𝜃 is the angle formed at atom (i) by neighbors (j) and (k). This quantity is\n𝑗𝑖𝑘\ninvariant under rotation (𝑙 = 0) and corresponds to projecting the 𝑙 = 1 vectors onto\nthe scalar subspace via Clebsch-Gordan coupling in spherical harmonic space:\n∑1 𝑌𝑚(𝑟 ̂ )𝑌𝑚∗(𝑟 ̂ ) ∝ cos𝜃 . (10)\n𝑚=−1 1 𝑖𝑗 1 𝑖𝑘 𝑗𝑖𝑘\nThus, second-order invariants efficiently encode angular correlations while\nremaining rotationally consistent5.\nThird-order composite tensors further capture higher-order angular correlations:\n3 3\n(11)\n|𝑢⃗⃗⃗ |3 = (𝐴+𝐵)3/2 ≈ 𝐴3/2 + 3 𝐴1/2𝐵+ 3 𝐴−1/2𝐵2,\n𝑖 2 8\n2\nwith 𝐴 = ∑ |𝑢⃗⃗⃗⃗⃗ | and 𝐵 = ∑ 𝑢⃗⃗⃗⃗⃗ ⋅𝑢⃗⃗⃗⃗⃗ . This higher-order invariant remains in the\n𝑗 𝑖𝑗 𝑗≠𝑘 𝑖𝑗 𝑖𝑘\n$l=0$ subspace and can be interpreted as a projection of tensor products of\n$l=1$ vectors onto a scalar invariant, capturing multi-body angular correlations akin to\ncontributions from higher-order spherical harmonics ( 𝑙 > 1 ) without explicitly\nevaluating the spherical harmonic basis 4, 6.\nFour-body dihedral descriptors are defined as:\n(12)\n∑ (𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗⃗⃗ )⋅(𝑢⃗⃗⃗⃗⃗ ×𝑢⃗⃗⃗⃗ ) ∝ cos𝜑,\n𝑘,𝑙∈𝒩(𝑖,𝑗) 𝑖𝑘 𝑖𝑗 𝑖𝑗 𝑗𝑙\nwhere φ is the dihedral angle. This descriptor is invariant under rotation (𝑙 = 0) and\ncaptures the torsional geometry of the local environment. Formally, it can be\nrepresented in spherical harmonic space via Clebsch-Gordan coefficients:\n∑ ⟨1𝑚 ,1𝑚 |00⟩𝑌\n𝑚1(𝑟⃗⃗⃗⃗ ̂\n×𝑟⃗⃗⃗ )𝑌\n𝑚2(⃗𝑟⃗⃗ ̂\n×𝑟⃗⃗⃗ ) ∝ cos𝜑.\n(13)\n𝑚1,𝑚2 1 2 1 𝑖𝑘 𝑖𝑗 1 𝑗𝑙 𝑖𝑗\nThe plane defined by 𝑟⃗⃗⃗ is essential and cannot be omitted, as it determines the\n𝑖𝑗\norientation of the cross products. This construction provides a computationally efficient\napproximation of higher-order spherical harmonic interactions while maintaining\nrotational invariance4, 7.\nIn summary, LiTEN’s descriptors adhere rigorously to SO(3) representation\nprinciples. First-order vectors (𝑙 = 1) capture directional information, inner and cross\nproducts yield scalar invariants (𝑙 = 0 ), and multi-body higher-order descriptors\napproximate projections of higher-order tensor products, analogous to spherical\nharmonic convolutions. This design ensures that the network is both physically\nconsistent under rotations and expressive enough to distinguish local conformations\nwith different symmetries or directional biases efficiently.\nWe have also incorporated the above content into the Methods section of the\nSupplementary Information (SI), under “Theoretical Connection Between the LiTEN\nFramework and Spherical Harmonic Convolutions.”\nSuggestions 1. I recommend that the authors explore the connection between the\ndesigned functional forms and the higher-order and coupling terms in class-II classical\nforce fields.\nResponse: We thank the reviewer for this insightful suggestion. Indeed, the functional\nforms employed in LiTEN bear conceptual similarities to the higher-order and coupling\nterms commonly found in class-II classical force fields. In classical force fields, class-\nII terms explicitly account for bond–angle–angle, angle–angle, and torsion–bond\ncouplings, often through quadratic or cubic expansions of internal coordinates, to better\ncapture anharmonicity and cooperative effects in molecular geometries.\nIn LiTEN, the multi-body descriptors—including the three-body quadratic and\ncubic tensors (|𝑢⃗⃗⃗ |2,|𝑢⃗⃗⃗ |3 ) and the four-body dihedral terms—serve a conceptually\n𝑖 𝑖\nanalogous role. The quadratic term (|𝑢⃗⃗⃗ |2 = ∑ 𝑢⃗⃗⃗⃗⃗ ⋅𝑢⃗⃗⃗⃗⃗ ) encodes angular correlations\n𝑖 𝑗,𝑘 𝑖𝑗 𝑖𝑘\nbetween pairs of neighbors, reminiscent of angle–angle coupling in class-II potentials.\nThe cubic expansion (|𝑢⃗⃗⃗ |3 ≈ 𝐴3/2 + 3 𝐴1/2𝐵+ 3 𝐴−1/2𝐵2) naturally incorporates self-\n𝑖\n2 8\nand cross-interaction terms, akin to higher-order couplings that account for\nanharmonicity in classical force fields. Finally, the four-body dihedral descriptors\ncapture torsional interactions and their dependence on local geometries, analogous to\ntorsion–bond and torsion–angle couplings.\nImportantly, while class-II force fields rely on predefined analytic functional forms\nand parameters for specific coupling terms, LiTEN learns these interactions directly\nfrom data through neural network embeddings. This allows the model to capture both\nstandard couplings and more subtle, nontrivial multi-body correlations without\nmanually specifying functional forms. The resulting descriptors thus bridge the gap\nbetween classical coupling concepts and learned many-body representations, providing\na data-driven, physically informed analog of class-II multi-body and cross terms.\nWe have added a discussion in the revised manuscript to highlight this connection\nand emphasize that LiTEN’s multi-body descriptors can be viewed as a generalized,\nlearnable extension of the higher-order couplings in class-II classical force fields.\nWe have added the relevant content to the Methods section of the Supplementary\nInformation (SI), under the subsection titled “Theoretical Connection Between the\nLiTEN Framework and Spherical Harmonic Convolutions.”\n“From a physical perspective, these descriptors serve a role analogous to class-II coupling terms in\ntraditional force fields. The quadratic invariant (|𝑢⃗⃗⃗ |2) parallels angle–angle couplings, while the\n𝑖\ncubic term (|𝑢⃗⃗⃗ |3) introduces anharmonic cross terms, similar to bond–angle–angle couplings. The\n𝑖\nfour-body dihedral descriptor naturally captures torsion–angle and torsion–torsion interactions.\nHowever, in contrast to classical force fields that rely on predefined analytic expressions, LiTEN\nlearns these couplings directly from data, providing a generalized, learnable extension of class-II\nmulti-body terms.”\nSuggestion 2. The equations 15-20 are written somewhat unclearly with confusing\nnotations. The readability would be greatly improved if clearer mathematical expressions\nare used.\nResponse: Thank you for your suggestion. Equations 15-20 indeed appeared somewhat\nconfusing, so we have reorganized the formulas. The revisions have been incorporated\ninto the updated manuscript, particularly in the 'Tensorized Quadrangle Attention'\nand 'Scalar-Vector Fusion' sections of the Methods:\n“To incorporate node-level information and further enrich the representation, we apply a\nlearned attention mechanism. Specifically, attention scores are computed as:\n𝑎𝑙 = 𝑆𝑖𝐿𝑈 (ℎ𝑙 +𝐸𝑙 +ℎ𝑙)⊙ 𝛼, (15)\n𝑖𝑗 𝑖 𝑖𝑗 𝑗\nwhere ℎ𝑙 ∈ 𝑅𝐶 and ℎ𝑙 ∈ 𝑅𝐶 are the scalar features of atoms 𝑖 and 𝑗, and 𝛼 ∈ 𝑅𝐶 is\n𝑖 𝑗\na trainable parameter vector that modulates attention weights. This attention 𝑎𝑙 ∈ 𝑅𝐶 is then\n𝑖𝑗\ngated by a continuous cutoff function based on the interatomic distance:\n𝑎𝑙 = 𝑎𝑙 ⊙𝑅𝑎𝑑𝑖𝑎𝑙_𝑓𝑢𝑛𝑐𝑡𝑖𝑜𝑛(𝑑 ). (16)\n𝑖𝑗 𝑖𝑗 𝑖𝑗\nThis hybrid design enables adaptive feature weighting that respects both learned chemical\nrelevance and geometric locality. The updated node and vector features are then computed as:\n𝑠𝑐𝑎 𝑙 = ℎ𝑙 ⊙𝐸𝑙 ⊙𝑎𝑙 ,𝑣𝑒𝑐 𝑙 = 𝑢𝑙 ⊙ℎ𝑙 +𝑢𝑙 ⊙ℎ𝑙. (17)\n𝑚𝑖𝑗 𝑗 𝑖𝑗 𝑖𝑗 𝑚𝑖𝑗 𝑗 𝑗 𝑖𝑗 𝑗\nThe expression 𝑠𝑐𝑎 𝑙 ∈ 𝑅𝐶 represents scalar message passing, while 𝑣𝑒𝑐 𝑙 ∈\n𝑚𝑖𝑗 𝑚𝑖𝑗\n𝑅𝐶represents vector message passing. Finally, both scalar and vector features are aggregated across\nneighbors using scatter-based summation and residual connections:\nℎ𝑙 =scatter_sum (𝑠𝑐𝑎 𝑙 )+ℎ𝑙, 𝑢𝑙 =scatter_sum(𝑣𝑒𝑐 𝑙 )+𝑢𝑙. (18)\n𝑚𝑖𝑗 𝑚𝑖𝑗\nThe expression ℎ𝑙 ∈ 𝑅𝐶 represents the scalar features of the node at the 𝑙 -th layer ,\nrespectively. Similarly, 𝑢𝑙 ∈ 𝑅𝐶×𝟛 represents the vector features of the node at the 𝑙 -th layer.\nThis construction results in a fully vectorized, geometry-aware message passing scheme that models\ncomplex conformational dependencies in a scalable and differentiable manner. By capturing implicit\nfour-body interactions through local vector operations, Tensorized Quadrangle Attention enhances\nthe model's capacity to represent torsional energy landscapes and anisotropic molecular features, all\nwhile maintaining efficiency suitable for large-scale systems.”\n“This component serves to further integrate the vector and scalar node representations derived from\nthe Message Passing and TensorCrossFusion modules. The fusion mechanism is designed to leverage\nboth three-body and higher-order geometric interaction cues, whose theoretical underpinnings have\nbeen rigorously justified in the Model section. For each node 𝑖, the vector feature 𝑢𝑙 ∈ 𝑅𝐶×𝟛and\n𝑖\nscalar feature be ℎ𝑙 ∈ 𝑅𝐶 are processed to capture complex geometric correlations and reinforce\n𝑖\nthe bidirectional exchange of information between vectorial and scalar modalities. To this end, the\nvector representation is first projected into two separate subspaces:\n[ 𝑢 𝑙(1) , 𝑢 𝑙(2) ] = 𝐿𝑖𝑛𝑒𝑎𝑟 ( 𝑢𝑙 ). (19)\n𝑖 𝑖 𝑣𝑒𝑐 𝑖\nThe expressions 𝑢 𝑙(1) ,𝑢 𝑙(2) ∈ 𝑅𝐶×𝟛 represent the split node vectors in the latent space at\n𝑖 𝑖\nlayer 𝑙. These subspace projections are used to construct two geometric interaction terms: a vector\ntrilinear term, which captures directional alignment (akin to three-body angular interaction), and a\nvector multi-body term, which encodes magnitude-based modulation (analogous to multi-body\ngeometric scaling). These terms are formally defined as:\n3\n𝑙(1) 𝑙(2)\n𝑉𝑒𝑐𝑇𝑟𝑖 = ∑ 𝑢 ⊙ 𝑢 , (20)\n:,𝑗,∶ :,𝑗,∶\n𝑗=1\n3\n3\n2\n𝑉𝑒𝑐𝑀𝑢𝑡𝑖 = (√∑( 𝑢 𝑙(2) ) +𝜀) , 𝜀 = 10−8, (21)\n:,𝑗,:\n𝑗=1\nwhich 𝑉𝑒𝑐𝑇𝑟𝑖, 𝑉𝑒𝑐𝑀𝑢𝑡𝑖 ∈ 𝑅𝐶. These vector-derived signals are used to modulate scalar\nfeatures. The scalar representation is projected as:\n[ℎ𝑙(1),ℎ𝑙(2),ℎ𝑙(3)] = 𝐿𝑖𝑛𝑒𝑎𝑟 (ℎ𝑙). (22)\n𝑠𝑐𝑎\nThen, we compute the update residues:\nScalar update residue:\nℎ𝑙+1 = (𝑉𝑒𝑐𝑇𝑟𝑖 +𝑉𝑒𝑐𝑀𝑢𝑡𝑖)⊙ℎ𝑙(1) +ℎ𝑙(2). (23)\nVector update residue:\n𝑢𝑙+1 = 𝑢𝑙(1) ⋅ℎ𝑙(3). (24)\nwhere the scalar ℎ 𝑙(3) ∈ 𝑅𝐶is broadcast across the spatial dimension of 𝑢 𝑙(1) .\n𝑖 𝑖\nThese update residues are added back to the original features in subsequent stages, enabling\ntightly coupled scalar-vector representations that are sensitive to both directional geometry and\nfeature semantics.”\nReferences\n1. Wang X, Wu X, Brooks BR, Wang J. Accurate Free Energy Calculation via Multiscale\nSimulations Driven by Hybrid Machine Learning and Molecular Mechanics Potentials.\nJournal of Chemical Theory and Computation 21, 6979-6987 (2025).\n2. Xiao J, Chen Y, Zhang L, Wang H, Zhu T. A machine learning-based high-precision\ndensity functional method for drug-like molecules. Artificial Intelligence Chemistry 2,\n100037 (2024).\n3. Thomas N, et al. Tensor field networks: Rotation-and translation-equivariant neural\nnetworks for 3d point clouds. arXiv preprint arXiv:180208219, (2018).\n4. Batzner S, et al. E(3)-equivariant graph neural networks for data-efficient and accurate\ninteratomic potentials. Nature Communications 13, 2453 (2022).\n5. Schütt K, Unke O, Gastegger M. Equivariant message passing for the prediction of\ntensorial properties and molecular spectra. In: Proceedings of the 38th International\nConference on Machine Learning (eds Marina M, Tong Z). PMLR (2021).\n6. Anderson B, Hy TS, Kondor R. Cormorant: Covariant molecular neural networks.\nAdvances in neural information processing systems 32, (2019).\n7. Kondor R, Son HT, Pan H, Anderson B, Trivedi S. Covariant compositional networks\nfor learning graphs. arXiv preprint arXiv:180102144, (2018)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have adequately addressed my previous comments and have supplemented the manuscript with clearer\ndescriptions and additional results, including simulations of bulk water and explicitly solvated ala5-mer. I don't quite get the\npoint \"Furthermore, we extend this framework to many-body interactions\" and Eq. 6 in the author's response to my first\ncomment, but that's fine. The manuscript should be accepted for publication and will make a valuable contribution to the\nMLP field.", "reviewer_2": "The authors have carefully responded to my previous suggestions and I am satisfied with their responses.\nI was able to access the code and corresponding data and easily follow the provided instructions to run one of the training\nscripts."}, "author_response": ""}]

Decision letter