Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport

Zheyong Fan; Zezhu Zeng; Cunzhi Zhang; Yanzhou Wang; Haikuan Dong; Yue Chen; Tapio Ala-Nissila

Source-linked AI summary

Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport

Zheyong Fan, Zezhu Zeng, Cunzhi Zhang, Yanzhou Wang, Haikuan Dong, Yue Chen, Tapio Ala-Nissila

arXiv:2107.08119v2Jul 16, 2021physics.comp-ph

TL;DR

Existing neural-network potential training can be difficult because gradient descent may become trapped in local minima, while heat transport in strongly anharmonic or disordered materials is challenging to treat accurately. The paper develops NEP, combining evolutionary neural-network training, polynomial atomic-environment descriptors, and GPU implementation in GPUMD. NEP achieves accuracy comparable to other ML potentials with high computational efficiency and provides per-atom heat currents for heat-transport simulations.

Problem
Neural-network potentials have many fitting parameters, and gradient-descent training can become trapped in local minima; heat transport in strongly anharmonic or spatially disordered materials is also difficult to treat accurately.
Method
NEP uses a separable natural evolution strategy to train neural-network potentials with a Chebyshev- and Legendre-based atomic-environment descriptor, implemented on GPUs in GPUMD.
Results
NEP achieves accuracy comparable to GAP-SOAP, MTP, and DP while providing high computational efficiency, including about 0.1 microsecond per atom per MD step on one Nvidia V100 GPU.
Takeaways & Limitations
Per-atom heat currents combined with HNEMD and spectral decomposition enable efficient and accurate heat-transport simulations from quantum-mechanical training data.

Abstract

from arXiv · show

We develop a neuroevolution-potential (NEP) framework for generating neural network based machine-learning potentials. They are trained using an evolutionary strategy for performing large-scale molecular dynamics (MD) simulations. A descriptor of the atomic environment is constructed based on Chebyshev and Legendre polynomials. The method is implemented in graphic processing units within the open-source GPUMD package, which can attain a computational speed over $10^7$ atom-step per second using one Nvidia Tesla V100. Furthermore, per-atom heat current is available in NEP, which paves the way for efficient and accurate MD simulations of heat transport in materials with strong phonon anharmonicity or spatial disorder, which usually cannot be accurately treated either with traditional empirical potentials or with perturbative methods.

I. INTRODUCTION

ML potentials can approach quantum-mechanical accuracy at lower computational cost, but training neural networks is difficult because gradient descent may converge to suboptimal local minima. NEP addresses this with evolutionary training and GPU acceleration for large-scale MD and heat-transport simulations.

ML potentials can achieve accuracy comparable to quantum-mechanical training data while requiring only a small fraction of the computation time.
Neural-network potentials contain many fitting parameters, making parameter optimization difficult and gradient descent vulnerable to local minima.
Neuroevolution uses global-search evolutionary algorithms to optimize neural networks, including scalable separable natural evolution strategies.
NEP combines separable natural evolution strategy training with an efficient GPU implementation in the open-source GPUMD package.
The paper evaluates NEP against GAP, MTP, and DP implementations and focuses particularly on heat-transport simulations.
Many-body potentials require explicit per-atom energy, force, virial, and heat-current formulations for efficient implementation and simulation.

B. Virial stress, heat current and thermal conductivity

The framework derives per-atom virial and heat-current quantities from partial forces, enabling stress-controlled MD and detailed thermal-transport calculations. These quantities support conductivity and spectral-conductivity evaluations through HNEMD.

The per-atom virial is generally nonsymmetric, so heat-transport applications must use its full 3 × 3 tensor.
The reported heat-current expression has been extensively documented, whereas the LAMMPS implementation has been shown erroneous for general many-body potentials.
Total heat current is obtained by summing the per-atom heat currents, and HNEMD uses it to compute the lattice thermal-conductivity tensor.
HNEMD applies a driving-force parameter that induces an external force on each atom, while spectral conductivity resolves conductivity by angular frequency.
Explicit partial forces enable convenient calculation of per-atom virials and heat currents for constant-stress ensembles and heat-transport simulations.

1. Single-component systems

NEP represents atomic environments with radial and angular descriptor components built from Chebyshev and Legendre polynomials. Cutoffs, expansion orders, and weighting choices control the descriptor’s treatment of local environments and interaction ranges.

NEP’s descriptor contains radial components and angular components indexed by radial order n and angular order l.
Angular descriptor terms use Legendre polynomials of bond angles, while radial functions use first-kind Chebyshev polynomials of a transformed distance variable.
The transformed Chebyshev variable spans −1 to 1, and positive-definite radial functions were found empirically beneficial.
The cutoff function vanishes beyond rc, and the first four lowest-order radial functions are illustrated in Figure 1.
Radial and angular expansions use separate maximum radial orders, nR_max and nA_max, and can have different cutoff radii.
Radial components can represent relatively long-ranged interactions, whereas angular components mainly account for intermediate-ranged interactions.
A modified radial weighting produces a descriptor similar to SOAP, although the radial expansion schemes remain different.

2. Multicomponent systems

NEP incorporates atomic-species information into permutation-symmetric descriptors for multicomponent systems, including the central atom species to distinguish otherwise identical environments.

2. Multicomponent systems: Multicomponent NEP descriptors incorporate atomic species through weighting factors applied to radial and angular terms.The described approach uses factors such as Z_j and Z_jZ_k, with corresponding modification of the cutoff function.
2. Multicomponent systems: The modified descriptor reduces the relative influence of neighboring atoms with very different atomic numbers.The factor Z_j is changed so relative weights become less important when atomic numbers differ substantially.
2. Multicomponent systems: Including the central atom factor √Z_i distinguishes configurations with identical environments but different central atom species.Earlier weighting approaches could not distinguish those configurations.
2. Multicomponent systems: In NEP, the site energy is modeled as a function of descriptor components using a feedforward neural network.The descriptor components provide the variables of the scalar site-energy function.

1. The neural network model

The NEP neural network maps atomic-environment descriptors to site energy through hidden-layer transformations, while the resulting network supports force-related derivatives through descriptor expressions.

1. The neural network model: The descriptor vector serves as the neural network input and the site energy as its output.The model may contain one or more hidden layers, with a single hidden layer used for presentation.
1. The neural network model: Descriptor components are indexed as ν = (nl), forming a vector with N_des components.The combined index simplifies presentation of the descriptor-to-hidden-layer coupling.
1. The neural network model: The hidden-layer state is obtained from the descriptor vector through linear and nonlinear transformations.The hidden layer contains N_neu neurons, and the nonlinear activation is specified as hyperbolic tangent.
1. The neural network model: The site energy is computed as a linear combination of hidden-layer states using output weights and a bias.The cited passages identify the connection weights and output-layer bias in this calculation.
1. The neural network model: Descriptor derivatives with respect to interatomic coordinates provide the radial and angular components needed for force calculations.The force-related derivative terms are calculated from the descriptor expressions.

E. Training the machine learning potential

NEP training optimizes neural-network weights and biases by minimizing a loss that combines energy, force, virial, and regularization terms.

E. Training the machine learning potential: Training determines neural-network weights and biases by minimizing errors against quantum-mechanical energy, force, and virial data.The parameter vector z contains N_par fitting parameters.
E. Training the machine learning potential: The training dataset typically contains a few thousand structures in the examples studied.The passage states this number in the context of the loss-function formulation.
E. Training the machine learning potential: Training is expressed as a real-valued optimization problem seeking an optimal parameter set z* that minimizes the loss.The parameter vector and its optimum are introduced as the variables of the optimization.
E. Training the machine learning potential: The total loss is a weighted sum of energy, force, virial, L1, and L2 loss terms.Energy, force, and virial losses are defined using RMSE-based terms, while L1 and L2 provide regularization.
E. Training the machine learning potential: L1 and L2 regularization encourage smaller weight magnitudes and can help prevent over-fitting.The cited formulation applies both regularization types to the neural network.
E. Training the machine learning potential: For energy, force, and virial expressed in the stated units, λ_e = λ_f = λ_v = 1 is reported as a good choice.The regularization weights λ_1 and λ_2 must be tuned to balance over-fitting and under-fitting.

2. Separable natural evolution strategy as the training algorithm

NEP uses separable natural evolution strategy to train neural-network parameters by evolving a search distribution over candidate solutions across generations.

2. Separable natural evolution strategy as the training algorithm: Separable natural evolution strategy updates a search distribution by following the natural gradient of the loss function.It is presented as a derivative-free black-box optimizer for real-valued evolutionary optimization.
2. Separable natural evolution strategy as the training algorithm: The algorithm initializes a search distribution with mean m and standard deviation s in a solution space of dimension N_par.The initial mean components may be random values between −1/2 and 1/2, while standard deviations may be set to 0.1.
2. Separable natural evolution strategy as the training algorithm: Each generation samples a population of candidate parameter vectors from the current search distribution.The population uses Gaussian random vectors and has size N_pop.
2. Separable natural evolution strategy as the training algorithm: Candidate solutions are evaluated with the loss functions and sorted from small to large loss.This ranking supplies the basis for subsequent population updates.
2. Separable natural evolution strategy as the training algorithm: Rank-based utility values evolve the population toward better individual solutions.The utilities u_k are used in the natural-gradient update.
2. Separable natural evolution strategy as the training algorithm: The search distribution's mean and standard deviation are updated elementwise using learning rates η_m = 1 and η_s = (3 + log N_par) / 5.The mean and standard deviation control subsequent sampling of candidate solutions.

F. Computer implementation

The NEP framework is implemented in the open-source GPUMD package for training potentials, running molecular dynamics, and calculating selected phonon properties on GPUs.

F. Computer implementation: NEP-related calculations are implemented in the open-source gpumd package using CUDA, with separate executables for training, molecular dynamics, and phonon-property calculations.The nep, gpumd, and phonon executables support these respective tasks.
F. Computer implementation: The GPU implementation derives explicit per-atom quantities and maps one atom to one CUDA thread for efficient parallelization.This implementation closely follows the GPU implementation of the Tersoff potential.
F. Computer implementation: Performance is evaluated against GAP-SOAP, MTP, and DP implementations using case studies designed to assess machine-learning potentials across material phases.The comparison uses quip, mlip, and DeePMD-kit packages.
F. Computer implementation: DFT calculations provide training data for bulk PbTe and 2D silicene, while the general-purpose silicon data are taken from a published reference.The inputs and outputs of the nep executable are openly available in Zenodo.
F. Computer implementation: The general-purpose silicon dataset contains 2475 structures spanning isolated, solid, bonded, liquid, amorphous, surface, vacancy, and other defective configurations.Every structure has energy data, but virial data are not available for all structures.
F. Computer implementation: The 2D silicene dataset contains 914 rectangular 60-atom cells covering temperatures from 100 K to 900 K and biaxial strains from −1% to 1%.These cells were obtained using active learning implemented in the mlip package.
F. Computer implementation: The silicene electronic-structure calculations use the PBE functional, an optimized norm-conserving Vanderbilt pseudopotential, a 40 Ry cutoff, and a 3×3×1 k-point mesh.The electronic self-consistent loop used a 10^-10 Ry threshold.

3. Training dataset for bulk PbTe

The bulk PbTe training set combines DFT-MD and displaced low-temperature configurations, while NEP training and evaluation use tunable descriptors, regularization, and comparisons with other potentials.

3. Training dataset for bulk PbTe: The bulk PbTe dataset contains 325 triclinic 250-atom cells, including 305 DFT-MD structures from 300 K to 900 K and 20 randomly displaced cells.The displaced cells sample the low-temperature phase space of PbTe.
3. Training dataset for bulk PbTe: NEP hyperparameters include radial and angular cutoffs, Chebyshev and Legendre expansion orders, hidden-layer size, regularization parameters, population size, and generation count.Table I defines the roles of these tunable quantities.
3. Training dataset for bulk PbTe: For covalent silicon and silicene, the radial and angular descriptor cutoffs are not differentiated, whereas PbTe benefits from a larger radial and smaller angular cutoff.This choice reduces computational cost at a given target accuracy for PbTe.
3. Training dataset for bulk PbTe: A single hidden layer with 40 or 50 neurons is sufficient for all materials tested, unlike the deeper networks typically used with simpler descriptors in DP.The descriptor dimensions are 80, 65, and 41 for general silicon, silicene, and bulk PbTe, respectively.
3. Training dataset for bulk PbTe: With increasing generations, energy, force, and virial RMSEs decrease and converge, while L1 and L2 regularization losses first rise and then fall.The regularization prevents uncontrolled growth of neural-network weights and biases that can lead to over-fitting.
3. Training dataset for bulk PbTe: Independent training runs produce comparable results with very small variation, suggesting stable optimization by the natural evolution strategy.The passage frames this as evidence that the method finds good minima of the loss function.
3. Training dataset for bulk PbTe: Accuracy comparisons report energy and virial RMSEs in meV/atom and force RMSEs in meV/Å for NEP, GAP-SOAP, MTP, and DP.The comparisons are summarized in Table II and Figure 3.

D. Evaluation of the computational efficiency

NEP combines GPU execution with competitive accuracy and substantially higher computational efficiency than several ML potentials, supporting large-scale MD simulations and heat-transport studies.

Benchmark setup: Computational speed is measured in atom-step per second using fixed-size MD simulations and specified CPU or GPU resources.The benchmark uses 8 × 10^3 atoms for bulk silicon and PbTe, 6 × 10^3 atoms for silicene, and 100 MD steps.
Accuracy–speed trade-off: GAP-SOAP achieves very high accuracy for general silicon but has very low computational speed, while MTP is about one order of magnitude faster at comparable accuracy.The GAP-SOAP configuration uses nmax = 10, lmax = 12, and 9 × 10^3 basis functions; reducing these parameters increases speed but reduces accuracy.
NEP performance: More than 10^7 atom-step per second is attainable with NEP on one Nvidia Tesla V100 GPU.This speed is about one order of magnitude below efficient empirical potentials such as Tersoff in gpumd.
Comparison with other ML potentials: NEP is more than one order of magnitude faster than DP when both use the same V100 GPU.The DP comparison uses the recent model-compression technique to accelerate calculations.
Heat-transport application: NEP yields a silicene lattice thermal conductivity of 33.7 ± 0.6 W/mK, consistent with GAP-SOAP and within the BTE-DFT range.The GAP-SOAP value is 32.4 ± 2.9 W/mK, while both Stillinger-Weber values are well below the BTE-DFT range.
Physical validation: NEP phonon dispersions agree well with DFT, whereas both Stillinger-Weber versions significantly underestimate phonon frequencies near Γ.The frequency underestimation implies lower phonon group velocities around Γ for the Stillinger-Weber potentials.

B. Thermal transport in bulk PbTe

NEP-based simulations reproduce reported heat-transport behavior across crystalline and disordered silicon and bulk PbTe, while highlighting the importance of anharmonic scattering and quantum corrections.

Bulk PbTe: NEP predictions for bulk PbTe thermal conductivity from 300 to 700 K agree well with experimental data.Comparable results are also obtained by BTE calculations including both three-phonon and four-phonon scatterings.
Bulk PbTe: BTE calculations for PbTe considering only three-phonon scattering produce lattice thermal conductivities about two times larger than calculations including four-phonon scattering.The comparison indicates the importance of four-phonon scattering in PbTe.
Methodological scope: NEP heat-transport simulations are presented as promising for materials with strong phonon anharmonicity because MD fully accounts for phonon anharmonicity.The paper contrasts this scope with perturbative BTE treatments in strongly anharmonic materials.
Amorphous silicon: NEP reproduces the radial and angular distribution functions of liquid and amorphous silicon in agreement with DFT and GAP results.These structural comparisons support the reasonable structural properties of the generated amorphous-silicon samples.
Amorphous silicon: In amorphous silicon, low-frequency phonons mainly contribute to thermal conductivity despite high-frequency phonons contributing substantially to the phonon density of states.The result is described as similar to the corresponding behavior in crystalline 2D silicene.
Amorphous silicon: The amorphous-silicon thermal conductivity is 1.5 ± 0.1 W/mK classically and 1.4 ± 0.1 W/mK after quantum correction at 300 K.The quantum effect is minor at 300 K but can be very strong at low temperatures.

Appendix A: The Chebyshev polynomials used in NEP

The supplied appendix passages list Chebyshev-based MTP and GAP configurations and describe a silicene regularization test, but do not provide the NEP polynomial formulations themselves.

Polynomial definitions: The appendix states that Chebyshev polynomials of the first kind use initial values T0(x) = 1 and T1(x) = x with a recurrence relation.The supplied passage does not include the recurrence expression beyond identifying its existence for n ≥ 2.
Polynomial definitions: The derivative of Tn(x) is related to Chebyshev polynomials of the second kind Un−1(x) for n > 0.The supplied excerpt does not include the explicit derivative equation.
Polynomial definitions: The Chebyshev polynomials of the second kind use initial values U0(x) = 1 and U1(x) = 2x with a recurrence relation for n ≥ 2.Only the initial values and existence of the recurrence are stated in the supplied passage.
Appendix configurations: The appendix provides input commands for training MTP and GAP potentials, including Chebyshev radial-basis settings and SOAP configurations.Listed parameters include radial-basis sizes, cutoff distances, sparse basis counts, and energy or force fitting settings.
Regularization test: A silicene regularization test compares training on 714 structures with testing on 200 structures that include pressures and temperatures outside the training set.The test therefore includes both interpolation and extrapolation.
Regularization test: Increasing regularization from λ1 = λ2 = 0 to 0.05 and 0.5 reduces mean absolute neural-network parameter values from 12 to 0.36 and 0.12.The passage labels 0.05 as proper regularization and 0.5 as strong regularization.