1. 搜索关键词

• GGA, PBE, CP2K
• CP2K, glass

2. 论文计算方法摘要

1. The mobility and solvation structure of a hydroxyl radical in a water nanodroplet: a Born–Oppenheimer molecular dynamics study

• https://pubs.rsc.org/en/content/articlehtml/2021/cp/d1cp01830b

Computational methods

DFT-based Born–Oppenheimer molecular dynamics (BOMD) simulations were carried out using the Quickstep module of the CP2K package28 in which a Gaussian-typed basis is used to describe the wave functions while an auxiliary plane-wave basis is used to describe the electron density.29 To account for the unpaired electron of OH*, spin-polarized DFT calculations were performed using the BLYP exchange–correlation functional.30,31 Short MD simulations with the B3LYP hybrid functional were performed to check the effect of self-interaction error in the BLYP functional. Grimme's dispersion correction DFT-D3 with a zero damping function32 was applied to account for weak dispersion interactions. A double zeta valence polarized (DZVP) Gaussian basis set33 was chosen with the Goedecker–Teter–Hutter (GTH) norm conserving pseudopotentials34,35 for core electrons. BOMD simulations were performed under the canonical (NVT) ensemble with a 1 fs time step. The temperature was set at 300 K via coupling to a Nose–Hoover36 thermostat chain, with a time constant of 0.1 ps.

The simulation system consists of one OH* radical and 191 water molecules, forming a nearly spherical shape with a ∼10 Å radius. This water droplet was placed in the center of a 40 × 40 × 40 Å3 cubic box, which allows the formation of an air–water interface and prevents it from interacting with its adjacent periodic images. Two different initial positions for OH* were adopted, where the distances between OH* and the center of mass of the water droplet (RCOM-OH*) are 10.3 Å (at the surface) and 1.5 Å (in the interior region), respectively (Fig. 1). Henceforth, the two resulting trajectories are labelled as Sint and Ssur, respectively. The total simulation time was 150 ps for both trajectories.

To identify the relative position where the radical was located with respect to the air–water interface, Voronoi analyses were performed using “Voro++” software.37 For a specific configuration, the coordinates of all the oxygen atoms (including the radical oxygen) were used to form a three-dimensional point set. After performing Voronoi tessellation in the 40 × 40 × 40 Å3 cubic box, positions of each oxygen atom are classified into three categories. Oxygen atoms whose Voronoi cell shares a boundary with the cubic box (with its Voronoi cell volume larger than 30 Å3) are assigned to the surface layer. Oxygen atoms with at least one surface-layer neighbour are assigned to the subsurface layer. The remaining oxygen atoms are assigned in the inner (interior) region of the water droplet.

To quantify the free energy barrier of transferring a hydrogen atom from a neighbouring water molecule to the hydroxyl radical, slow-growth free energy calculations38 were performed using the difference between the O–H and O*–H distances as the constraint variable R. Hereafter, O and H are used to denote the oxygen and hydrogen atoms from water molecules, while O* and H* are used to denote those from the hydroxyl radical, respectively. The average constraint force was measured over a 6.4 ps trajectory with a 0.00025 Å fs−1 growth rate, allowing the R to span from −0.8 Å to 0.8 Å. The free energy profile was obtained from a thermodynamic integration over the defined coordinate R.

2. Ab Initio Molecular Dynamics of CdSe Quantum-Dot-Doped Glasses

• https://pubs.acs.org/doi/full/10.1021/jacs.9b12073

2.1. Generation of a Glass Matrix Structure

The starting configuration for the glass matrix was generated by placing atoms randomly in a cubic simulation box. The total number of atoms for the glass was 540 (120 Na, 120 Si, and 300 O), with the simulation cell sizes (a = b = c = 19.383 Å, α = β = γ = 90°) kept constant throughout the simulation, giving a density consistent with experimental values (ρ = 2.492 g/cm3). (15) Hard constraints were imposed to avoid unphysically small interatomic distances. An initial classical molecular dynamics simulation was performed using a partial-charge rigid-ion pairwise potential developed by Pedone et al., (16) with the DL_POLY classic package. (24) The Coulomb interactions were calculated using the Ewald summation method (25) with a precision of 10–5 and a real-space cutoff for short-range interactions set to 7.6 Å. The Verlet algorithm was applied for the integration of the equations of motion with a time step of 1 fs. The glass structures were generated using a melt-quenching approach in the NVT ensemble at the target density from experimental data, using a Nosé–Hoover thermostat (26−28) with a relaxation time of 0.1 ps. The initial structure was heated up gradually in steps of 1000 K with a 60 ps MD run at each temperature from 300 to 6000 K. After equilibration of the liquid at 6000 K during 400 ps, the system was cooled gradually in steps of 500 K with a 60 ps MD run at each temperature from 6000 to 300 K. Another 200 ps NVT simulation was carried out at 300 K, together with a 200 ps NVE simulation in order to equilibrate the structure.

2.2. Generation of the CdSe Quantum-Dot-Doped Glass Structure

The experimental mechanism for the formation of CdSe QDs in glass is based on the phase decomposition of oversaturated solid solution, (9) and the time scale of the growth of these nanocrystals is far beyond the current computational method. Therefore, it is difficult to simulate the nucleation and growth stages, and here we directly introduced the CdSe QD into the glass matrix rather than trying to model directly the melt-quenching process. The Cd33Se33 QD (with a diameter of 13 Å) obtained in our previous work (29) was incorporated into the glass matrix, removing glass atoms so interatomic distances between the QD and the glass matrix were longer than 2.5 Å. The final composition was 39 Na2O–78 SiO2–33 CdSe (Figure 1a). The interatomic interactions in CdSe QD used a Lenard-Jones pairwise potential validated in the literature (22,23) and the Lorentz–Berthelot combining rules were used for interactions between CdSe QD and the glass matrix. The whole structure was equilibrated at 1000 K, first using 200 ps of NVT dynamics and then 200 ps of NVE dynamics, initially with atoms frozen from the CdSe QD. The structure was then cooled gradually in steps of 50 K with a 60 ps MD run at each temperature from 1000 to 500 K while keeping CdSe QD frozen. Subsequent cooling from 500 to 300 K took place in steps of 10 K, with all atoms allowed to move. Finally, the structures were further equilibrated at 300 K with 200 ps NVT dynamics and final 200 ps NVE dynamics (Figure 1b). This simulation methodology was followed after comparison with several other simulations in order to allow full relaxation of the system (especially the glass at high temperature) while maintaining the integrity of the quantum dot. We have explored several quenching methods of the whole structures, and the size of the glass matrix has been altered to test the representation of our current model, which can be found in the Supporting Information.

2.3. Ab Initio Molecular Dynamics of a CdSe Quantum Dot in a Glass Matrix

After the equilibration of systems using the classical MD simulations described above, we used the resulting configurations as a starting point for ab initio modeling at the density functional theory (DFT) level, using the Kohn–Sham formulation as implemented in the CP2K code. (30) Simulations were run at the generalized gradient approximation level, employing the PBE exchange–correlation functional. (31) The plane wave cutoff was set to 600 Ry. For Na, Cd, and Se, we used short-range molecularly optimized double-ζ single polarized basis sets (DZVP-MOLOPT-SR-GTH), while for O and Si we used a double-ζ single polarized basis set (DZVP-MOLOPT-GTH). (32,33) After an initial geometry optimization with DFT, the resulting relaxed structure (Figure 1c) was used as the initial structure for the AIMD simulations. The structure was quenched from 500 to 300 K in steps of 50 K, with a total 10 ps AIMD run at a time step of 2 fs. The production run was conducted in the NVT ensemble at 300 K for 10 ps.

3. Impeding 99Tc(IV) mobility in novel waste forms

• https://www.nature.com/articles/ncomms12067

Methods

Density functional theory (DFT) parameters

Spin-polarized DFT simulations were performed with periodic boundary conditions (3D PBC) as implemented in the CP2K package37. The Perdew, Burke and Ernzerhof (PBE) generalized gradient approximation was used for the exchange-correlation functional38. The core electrons were described by the norm-conserving pseudopotentials39, while the valence wave functions were expanded in terms of double-zeta quality basis sets optimized for condensed systems to minimize linear dependencies and superposition errors40. An additional auxiliary plane wave basis set with a 500-Ry cutoff was used to calculate the electrostatic terms. The GGA+U scheme was used to provide more accurate electronic structure for the localized d-orbitals. The Hubbard parameter (U–J) of 3.5 eV was taken for the Fe 3d states, which results in a work function of 5.32 eV, in good agreement with that obtained by Pentcheva et al.28 Owing to large supercell simulations, the Γ-point approximation was used for the Brillouin zone integration.

Computational models

To study Tc incorporation in magnetite with and without dopants, we used a 2 × 2 × 2 supercell in all simulations to minimize periodic images. Optimization of the bulk structure of magnetite had a cell parameter of 8.391 Å, which agrees well with experimental data (8.390 Å (ref. 41)). Using this optimized cell parameter, we constructed a magnetite(001) surface model terminated at an octahedral Fe sublattice, since it is known to be the most stable surface structure in magnetite. A more recent surface model was also considered42, but was found not to have significant impact on the present problem, see SI. Our model system consisted of a symmetric slab with seven octahedral and six tetrahedral Fe sublattices (384 atoms) with a vacuum region of 12.5 Å between slabs. To study Tc incorporation, one surface octahedral Fe was replaced with Tc, followed by structural optimization. We also optimized a structure with one octahedral Fe in the third layer replaced by Tc. For the doping studies, we substituted a surface Fe atom with Co, Ni or Zn (∼1 wt%) at a lattice position close to Tc. In all simulations, we fixed the atomic positions of the four bottom atomic layers.

AIMD simulations

AIMD simulations were performed with and without Tc at 25 °C and with the dopants Co/Ni/Zn at 600 °C, with the Nosé–Hoover thermostat for NVT ensemble and a time step of 1.0 fs. Each simulation was equilibrated for at least 20–28 ps, and the last 10–12 ps of the trajectories was used for the analysis. Owing to the big computational cost of high-temperature simulations, we chose lower range of vitrification temperatures (600 °C), while experiments were performed at somewhat higher temperatures (∼700 °C).

4. Influence of Glass Composition on the Luminescence Mechanisms of CdSe Quantum-Dot-Doped Glasses

• https://pubs.acs.org/doi/full/10.1021/acs.jpcc.1c04665

We generated configurations of pristine sodium silicate glasses and, from those pristine glasses, created CdSe quantum-dot-doped silicate glasses, using a combination of classical force-field-based molecular dynamics (MD) and ab initio molecular dynamics (AIMD). The protocol followed here is adapted from that developed and validated in our previous work. (13) The pristine QD was cut from the bulk structure and then relaxed to the lowest energy configuration. (12) The detailed methodology involved in the generation of glass matrix and CdSe quantum-dot-doped glasses is presented in the Supporting Information.

For each system, 500 configurations were selected at equal intervals from a 10 ps production run of ab initio molecular dynamics simulations. Among those, 50 configurations were picked to perform the calculation and analysis of the electronic structure of CdSe quantum-dot-doped glasses. We used the nonlocal PBE0-TC-LRC exchange-correlation functional with a cutoff radius of 2.0 Å. (15) From previous work, the electronic structure calculations based on the generalized gradient approximation (GGA) or local density approximation (LDA) exchange-correlation functionals severely underestimated the HOMO–LUMO gap and long-range Coulombic interactions. (16−18) The inclusion of the Hartree–Fock exchanges was found in the literature, to provide a more accurate description (18,19) of the band gap. Moreover, the computational cost of nonlocal functional calculations can be reduced using the auxiliary density matrix method (ADMM) (Table 1). (20)

5. Modelling the local atomic structure of molybdenum in nuclear waste glasses with ab initio molecular dynamics simulations

• https://pubs.rsc.org/en/content/articlehtml/2016/cp/c6cp03076a

Methods

Two systems with composition (SiO2)57.5–(B2O3)10–(Na2O)15–(CaO)15–(MoO3)2.5 and (SiO2)57.3–(B2O3)20–(Na2O)6.8–(Li2O)13.4–(MoO3)2.5 (thereafter referred as MCNB and MLNB respectively) were modelled using periodic boundary conditions. The molar composition of the MCNB glass is based on nuclear waste glass compositions that have been previously studied experimentally.12 In order to investigate the dependence of Mo local atomic structure and its stability within the borosilicate glass network on the chemical nature and composition of the glass host, we increased the molar composition of B2O3, decreased the concentration of Na2O and substituted the CaO with Li2O keeping its concentration at similar level, resulting in the MLNB composition. The concentration of MoO3 was the same in both glasses in order to avoid effects due to compositional dependence of the molybdenum.

The starting configurations were generated by placing atoms randomly in a cubic simulation box with imposed constraints to avoid un-physically small interatomic distances. The total number of atoms for the MCNB glass was 202 (38 Si, 12 B, 20 Na, 10 Ca, 2 Mo and 120 O), while for the MLNB glass it was 226 (38 Si, 26 B, 8 Na, 18 Li, 2 Mo and 134 O). The simulation cell sizes, 13.61 Å and 13.90 Å for MCNB and MLNB respectively, were chosen to give the appropriate density and kept constant throughout the simulation.

The density was calculated for the alkali alkaline-earth and mixed alkali borosilicate glasses before the addition of the MoO3. Classical MD simulations with the DL_POLY classic package,26 using a pairwise Lennard-Jones potential model27–29 and a melt-and-quench approach, were performed in order to generate the amorphous structures. The NPT ensemble (constant number of particles, pressure and temperature), using a Nosé–Hoover thermostat and barostat,30–32 was applied to calculate the density of each structure at 300 K and 0 bar. The accuracy of the calculated density was subsequently increased by running AIMD simulations at 300 K with the NPT ensemble and a CSVR thermostat and barostat.33 The calculated density for the calcium sodium borosilicate glass structure (2.59 g cm−3) is in very good agreement with the experimental value (2.62 g cm−3).34 There are no experimental density measurements available for the exact composition of the simulated mixed alkali borosilicate glass, however, the calculated density of this glass composition (2.49 g cm−3) is in good agreement with experimental densities (2.43–2.47 g cm−3) for lithium sodium borosilicate glasses of similar molar composition.35

The densities of the Mo-containing nuclear waste glasses modelled in this study were not obtainable experimentally, nor available via glass property modelling databases. We have therefore estimated the effect on the density of incorporating Mo to the two glass compositions performing a cell optimisation using variable cell AIMD with 1.0 fs timestep and a convergence threshold of 0.1 kbar for the components of the stress tensor.

Born–Oppenheimer AIMD simulations were carried out using the CP2K code.36 The electronic structure was treated through the Kohn–Sham formulation of density functional theory (DFT)37 using the generalised gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) exchange–correlation functionals.38 The Gaussian basis set employed for all the atomic species was a double-ζ basis set with polarisation functions (DVZP)39 in conjunction with the Goedecker–Teter–Hutter (GTH) pseudopotential.40 The plane-wave energy cutoff was 700 Ry and the MD timestep was 2.0 fs.

The glass structures were generated using a melt-and-quench approach. The canonical ensemble (constant number of particles, volume and temperature or NVT) was applied and the Nosé–Hoover thermostat chain,30–32 with a relaxation constant 0.1 ps, was chosen to control the temperature fluctuations. For each composition, the initial configuration was heated up at 2300 K with a 25 ps AIMD run to ensure that the system was melted and well equilibrated at this temperature. Despite a small drift in the total energy the recorded energy fluctuations were lower than 0.001%. The molten structure was subsequently cooled using a stepwise process, consisting of a series of nine NVT AIMD runs of 10 ps each, with target temperatures set to 2000 K, 1800 K, 1600 K, 1400 K, 1200 K, 1000 K, 800 K, 600 K and 300 K. At 300 K the structure was further equilibrated for 10 ps, followed by a final AIMD production run of 10 ps, to collect the structural data. This computational scheme corresponds to a total simulation time of 135 ps and a nominal cooling rate of around 20 K ps−1. Cooling rates of this order of magnitude have been used in previous simulation studies, using AIMD,21,23–25,41,42 in order to prepare accurate structural models of glasses that are in agreement with experimental results.

We note that, the fixed volume approach for the melt-and-quench process will generate glass structure with high final pressure in the cell. However, this residual pressure does not affect the quality of the calculations, as relaxation of the cell results in almost no difference in the final glass volume for both compositions (see Tables S1 and S2, ESI†).

6. Effect of alkaline oxides on aluminate slag structure by first principles calculation

• https://doi.org/10.1016/j.molliq.2023.123088

Simulation methods and experimental procedure

CP2K 24(https://www.sciencedirect.com/science/article/pii/S0167732223018949#b0120) /Quickstep 25(https://www.sciencedirect.com/science/article/pii/S0167732223018949#b0125) was used for first-principles molecular dynamics simulation in present study. The 3 s, 3p, 4 s electrons of Ca, 2 s, 2p electrons of O, 3 s, 3p electrons of Al, 2 s, 2p, 3 s electrons of Na, 1 s, 2 s electrons of Li, 5 s, 5p, 6 s electrons of Ba and 2 s, 2p, 3 s electrons of Mg were treated as valence. More simulation details had been described in previous work 26(https://www.sciencedirect.com/science/article/pii/S0167732223018949#b0130). In order to study the effect of alkaline oxides in aluminate slag, the atomic number and density of each sample was shown in Table 1. Thereinto, AC, NAC, LAC, BAC, and MAC are abbreviations for Al2O3-CaO, Na2O-Al2O3-CaO, Li2O-Al2O3-CaO, BaO-Al2O3-CaO, and MgO-Al2O3-CaO, respectively. As shown in Fig. 1, the simulation process includes four stages: mixing uniformity, cooling, computational density, and statistical analysis. At the end of the calculation, all frames during the final 6.15 ps data of Step 4 were statistically analyzed, and then calculate the average value of the results. In order to ensure the correct statistics of the atomic structure at the boundary, the atoms of the simulation results were expanded into a 333 supercell, and then the structure formed by the initial cell atoms and the supercell atoms was analyzed. Moreover, the data accounting for less than 1 % were ignored in the statistical process.

To further verify the accuracy of simulation results, the structural characteristics of molten slag were determined by water quenching method from 1873 K to obtain the vitreous body with an approximately retained high-temperature structural state. The water quenched slag was placed in an oven at 378 K and dried for 24 h. The dried slag was ground to a particle size of ∼ 45 µm using an agate bowl. The structural characteristics of the slag network were evaluated using Fourier transform infrared (FTIR) spectroscopy (Thermo Scientific Nicolet 6700) and X-ray photoelectron spectrometer (XPS, Thermo Scientific K-Alpha). For FTIR analysis, the vitreous powder (2 mg) was mixed with metal halide KBr (200 mg) and ground evenly. The mixed powder was held under pressure (20 ∼ 30 MPa), and pressed into thin slices using a tablet press to obtain samples for FTIR analysis. The O1s XPS detection was carried out with an Al Kα focused monochromator (1486.6 eV). To correct for the charging effects, the C1s peak with binding energy at 284.8 eV was used as the calibration standard for the experimental slags 27(https://www.sciencedirect.com/science/article/pii/S0167732223018949#b0135).

7. Atomic and electronic structures of an extremely fragile liquid

• https://www.nature.com/articles/ncomms6892

Methods

Density measurement

The density measurement of l-ZrO2 was performed with an aerodynamic levitator6,47,48. A small ZrO2 sample whose diameter was around 2 mm was set in a shallow nozzle where the sample was aerodynamically levitated. The levitated sample was then heated by a 100-W CO2 laser and a 500-W Nd:YAG laser. The temperature of the sample was measured by a single colour pyrometer. The weight of the recovered sample was measured. The temperature was calibrated using the given melting temperature (Tm=2,715 °C) in density measurements. The details of measurement can be found in the Supplementary Note 3 and typical image of the levitated specimen and the density data is shown in Supplementary Figs 6 and 7, respectively.

High-energy X-ray diffraction measurement

The high-energy X-ray diffraction experiments of l-ZrO2 were carried out at the BL04B2 and the BL08W beamlines49 at the SPring-8 synchrotron radiation facility, using the aerodynamic levitation technique6,47,48,49. The energy of the incident X-rays was 113 keV (BL04B2) and 116 keV (BL08W). The ZrO2 sample of 2-mm size was levitated by dry air and heated by a 100-W CO2 laser. The temperature of the sample was monitored by a two-colour pyrometer. As can be seen in Supplementary Fig. 8, the background of the instrument was successfully reduced due to adequate shielding of the detector and the optimization of beam stop. The measured X-ray diffraction data were corrected for polarization, absorption and background, and the contribution of Compton scattering was subtracted using standard analysis procedures49. The corrected data sets were normalized to give the Faber–Ziman16 total structure factor S(Q) and the total correlation function T(r) was obtained by a Fourier transformation of S(Q).

DF–MD simulation

The combined DF and MD simulations were performed with the CP2K programme package ( http://www.cp2k.org/)[50](https://www.nature.com/articles/ncomms6892#ref-CR50). The CP2K method uses two representations for the KS orbitals and electron density: localized Gaussian and plane wave basis sets. For the Gaussian-based (localized) expansion of the KS orbitals, we used a library of contracted molecularly optimized valence double-zeta plus polarization basis sets51 and the complementary plane wave basis set for electron density has a cutoff of 400 Ry. The generalized gradient approximation of Perdew, Burke and Ernzerhof52 (PBE) was adopted for the exchange-correlation energy functional and the valence electron–ion interaction was based on the norm-conserving and separable pseudopotentials of the form derived by Goedecker et al.53 We consider the following valence configurations: Zr (Cl4s24p25s25d2) and O (2s22p4). Periodic boundary conditions were used, with a single point (k=0) in the Brillouin zone. Effective charges of individual atoms were evaluated from electron density by integrating electronic charge inside the corresponding atomic volumes24. For reference, electronic structure of the high-temperature phase of c-ZrO2 (T>2,370 °C) is computed for a sample of 324 atoms.

The initial atomic configuration is created by a reverse Monte Carlo (RMC) simulation with high-energy X-ray diffraction data on 501 atoms in a cubic box of 18.98 Å (experimental density, 0.0733 atoms per Å3). The RMC++ programme code54 was used. The Born–Oppenheimer MD simulations were performed with a Nóse–Hoover thermostat55 and time steps of 2 fs (initialization) and 1 fs (data collecting). The system was simulated at 3,100 K (~2,800 °C) for a total of 30 ps, where the last 10 ps were used for data collection55. The corresponding mean-square displacement of atoms shows clearly a liquid behaviour (diffusion). The comparison of the partial pair correlation functions, gij(r), between the initial RMC configuration (start) and the DF–MD simulation is shown in Supplementary Fig. 9. The sharp shape of O–O gij(r) in the RMC model is artificially sharp due to small weighting factor for X-rays, while the shape of O–O gij(r) is reasonable in the DF–MD simulations. The system lost its memory of the initial (RMC) starting structure within a few picoseconds (Zr–O bond lifetimes ~185 fs).

8. First-principles study on microstructure of CaO-Al2O3-B2O3 slag

• https://www.sciencedirect.com/science/article/pii/S0167732222022772

Simulation methods and experimental procedure

The first-principles molecular dynamics simulation were carried out by CP2K 28(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0140) /Quickstep 29(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0145). The orbital transformation procedure 30(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0150) was used for the wave function optimization. The 2 s, 2p electrons of B, 2 s, 2p electrons of O, 3 s, 3p electrons of Al and 3 s, 3p, 4 s electrons of Ca were treated as valence, and the rest core electrons were represented by Goedecker-Teter-Hutter (GTH) pseudopotentials 31(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0155). The Gaussian basis set was double-ζ with one set of polarisation functions (DZVP) 32(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0160), and the plane wave cutoff was set to 500 Ry. Perdew-Burke-Ernzerhof (PBE) 33(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0165) functional was used to describe the exchange–correlation effects, and the dispersion correction was applied in all calculations with the DFT-D3 34(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0170) method. The condition of energy convergence was 1.0E-6 Hartree, and only one K-point at Γ (0, 0, 0) was used for Brillouin-zone sampling. A Nose-Hover thermostat was used to control the temperature in calculation.

The atomic number of each sample in the study system was shown in Table 1, and the density of the initial structure is calculated according to the empirical formula in reference 35(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0175). As shown in Fig. 1, the simulation process was divided into four stages. In the first stage, the canonical (NVT) ensemble retention system was adopted to simulate 6 ps at 4000 K, so that the atoms were mixed uniformly. The second stage used the micro-canonical (NVE) ensemble for reduction of system temperature to 1873 K 36(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0180), which took about 1.7 ps. In the third stage, the constant-pressure, constant-temperature (NPT) ensemble was applied to maintain the system to simulate 30 ps at 1873 K, and the densities of different components were obtained. At the last stage, the NVT ensemble was used to keep the system running 6.3 ps at 1837 K. To improve the calculation efficiency, the time step of the first two stages was set to 2 fs, while the last two stages was set to 1 fs. At the end of the calculation, the final 6.15 ps data of the fourth stage were statistically analyzed, and the data accounting for less than 1 % were ignored in the statistical process.

To further verify the accuracy of simulation results, the structural characteristics of molten slag were determined by water quenching method from 1873 K to obtain the vitreous body with an approximately retained high-temperature structural state. The water quenched slag was placed in an oven at 378 K and dried for 24 h, then ground to a particle size of ∼ 45 µm using an agate bowl. The structural characteristics of the slag network were evaluated using Fourier transform infrared spectroscopy (FTIR) and X-ray photoelectron spectroscopy (XPS). For FTIR analysis, the vitreous powder (2 mg) was mixed with metal halide KBr (200 mg) and ground evenly. The mixed powder was held under pressure (20 ∼ 30 MPa), and pressed into thin slices using a tablet press to obtain samples for Fourier transform infrared spectrometer (Thermo Scientific Nicolet 6700) analysis. The O1s XPS detection was carried out with an X-ray photoelectron spectrometer (Thermo Scientific K-Alpha) equipped with Al Kα focused monochromator (1486.6 eV). To correct the charging effects, the C1s peak with binding energy at 284.8 eV was used as the calibration standard for the experimental slags 37(https://www.sciencedirect.com/science/article/pii/S0167732222022772#b0185).

9. Partition model for trace elements between liquid metal and silicate melts involving the interfacial transition structure: An exploratory two-phase first-principles molecular dynamics study

• https://www.sciencedirect.com/science/article/pii/S0167732222015860#f0015

2. Computational method

2.1. Two-phase first-principles molecular dynamics simulation

Ab initio molecular dynamics simulations are time-consuming compared with classical molecular dynamics simulations. Thus, it is crucial to choose software that can perform efficient electronic structure calculations and molecular dynamics simulations. The widely used software that can perform first-principles molecular dynamics simulations include VASP, QUANTUM ESPRESSO, CASTEP, and CP2K. Among these software packages, CP2K 38(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0190) is excellent open-source software that takes into account both computing speed and accuracy, and it has been widely used in cutting-edge research in recent years. The QUICKSTEP 39(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0195) module implemented in CP2K is robust and efficient due to novel algorithms. The GPW method 40(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0200) adopted by CP2K is based on the hybrid Gaussian plane-wave method instead of traditional plane-wave basis sets, enabling accurate density functional calculations and efficient molecular dynamics simulations for condensed phases at a significantly reduced cost. In the present study, the two-phase first-principles molecular dynamics simulations were carried out in the framework of density functional theory (DFT) by using QUICKSTEP module as implemented in the CP2K 8.1 program package. A hybrid Gaussian plane-wave (GPW) method was employed that combines double-zeta valence plus polarization basis sets (DZVP) 41(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0205) with the generalized gradient approximation of the Perdew-Burke-Ernzerhof (PBE) exchange–correlation functional 42(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0210) and Goedecker-Teter-Hutter(GTH) norm-conserving pseudopotentials 43(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0215), 44(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0220), 45(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0225). A plane wave cutoff of 600 Ry for the Kohn-Sham orbitals was chosen as a compromise between accuracy and computational speed.

In this study, the local coordination structures of the partitioning element needed to be derived from two-phase ab initio molecular dynamics simulations before using the structure-based partition model to predict the distribution ratio. The initial compositions of the silicon phase (Sisingle bondB) and slag phase (CaOsingle bondSiO2) should be determined before performing the simulations, as shown in Fig. 1(c). The partitioning behavior of the impurity element B between the silicon melt and molten slag with different basicity at 1823 K was simulated, and the corresponding experimental partitioning data were available. In the initial stage of the simulation, the atoms of the two phases were randomly and uniformly placed in the starting simulation box, and the size of the simulation box was determined by the density of the silicon and silicate melt. As ab initio molecular dynamics simulation proceeded, the system eventually segregated into two dynamically equilibrated phases: the silicon melt and silicate melt. Finally, the distribution ratio was calculated by the newly established partition model based on the local coordination structures of the partitioning element and compared with the experimental data.

For each system, an AIMD run was started from the initial random configuration in the NVT ensemble at 4500 K for 10 ps of simulation time to ensure that the system was completely melted and equilibrated at this temperature, as shown in Fig. 1(a) and (b). Then, the molten silicon–silicate coexistence system was cooled using a stepwise process including a series of AIMD runs 46(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0230), 47(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0235), 48(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0240), 49(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0245), 50(https://www.sciencedirect.com/science/article/pii/S0167732222015860#b0250). Each model was run for 10 ps in NVT ensembles at each of the following temperatures: 4000 K, 3500 K, 3000 K, 2700 K, 2400 K, 2150 K, 2000 K, 1900 K and 1823 K. The production run was conducted in the NVT ensemble at 1823 K for 15 ps, and the last 10 ps of the trajectories was used for the time-averaged structural analysis. The time step of the AIMD simulation in this study was chosen as 0.5 fs. Periodic boundary conditions were applied to all the simulation systems. All first-principles molecular dynamics simulations are performed on the Beijing Super Cloud Computing Center.

10. Local structure and dynamics of ions in LiCl-GdCl3, KCl-GdCl3 and LiCl-GdCl3-Gd2O3 melts: Ab initio molecular dynamics simulations and Raman spectroscopy

• https://doi.org/10.1016/j.molliq.2023.121485

Simulation details

The melts of 50% GdCl3 – 50% LiCl, 49% GdCl3 – 49% LiCl – 2% Gd2O3, 50% GdCl3 – 50% KCl was simulated at a temperature of 1273 K by the method of ab initio molecular dynamics. For this purpose, three ensembles were used: 24Gd + 24Li + 96Cl, 26Gd + 24Li + 96Cl + 3O and 24Gd + 24K + 96Cl; hereinafter they are designated as GL, GLO and GK, respectively. The molecular dynamics cells were of a cubic shape with enabled periodic boundary conditions. The initial ionic distributions were obtained by running classical molecular dynamics with the Born-Mayer pair potential: E(rij) = qiqj/rij + A·exp(-rij/ρ). Here, rij is the interionic distance, qi, qj are the ionic charges, A and ρ are the adjustable parameters of short-range repulsion. The repulsive parameters of the potential were chosen to be intentionally incorrect. Namely, A = 5000 eV, ρ = 0.3 Å for all ion pairs. On the one hand, this allows avoiding strong non-optimal local energy minimum when switching to the ab initio method. On the other hand, some reasonable charge-ordered distribution of ions could be obtained. It is important to note that gadolinium oxide was completely dissociated at the starting point of ab initio simulations. The temperature was controlled by a Nozier–Hoover thermostat 17(https://www.sciencedirect.com/science/article/pii/S016773222300288X#b0085) with a time constant of 100 fs. The step of molecular dynamics was 5 fs.

The density functional theory with the PBE functional 18(https://www.sciencedirect.com/science/article/pii/S016773222300288X#b0090) was used to calculate the energy and forces. Ions were described by a combination of split-valence DZVP basis sets and Goedecker-Teter-Hutter pseudopotentials 19(https://www.sciencedirect.com/science/article/pii/S016773222300288X#b0095). In our opinion, double-zeta basis set represent a reasonable compromise between flexibility and computational speed, while compatible Goedecker-Teter-Hutter pseudopotentials are available for the majority of chemical elements. This allows studying properties of a variety of systems within the same approach. Therefore, the results are directly comparable. Namely, the trends in properties could be associated with the chemical nature, but not with the difference in computational approaches. The planewave cutoff was 500 Ry. The calculations were carried out using the CP2K code 20(https://www.sciencedirect.com/science/article/pii/S016773222300288X#b0100).

The density of the GL and GK ensembles was determined by preliminary simulation at a constant pressure of 1 bar and a temperature of 1273 K. The length of the simulations here was 13 ps. The cell lengths reduced down to their equilibrium values for the time of 2.8 and 1.5 ps for GL and GK ensembles, correspondingly. The rest of the trajectory allows satisfactory averaging of density. The equilibrium density was found to be 3.763 g/cm3 for GL and 3.36 g/cm3 for GK. Subsequently, these values ​​were used for modeling at constant volume. Taking into account that the content of gadolinium oxide is small, the density of the GLO ensemble was determined in an additive manner, using the crystallographic density of gadolinium oxide of 7.62 g/cm3 21(https://www.sciencedirect.com/science/article/pii/S016773222300288X#b0105). Thus, the density of the GLO ensemble was 3.84 g/cm3.

To study the structure and ion dynamics, simulations were carried out at a constant predetermined density for 50 ps. Due to preliminary modeling at constant pressure, the structure of the melts is well optimized, and the potential energy reaches an equilibrium value within several tens of molecular dynamic (MD) steps. This short starting interval was ignored for further analysis.

11. Liquid metal–organic frameworks

• https://www.nature.com/articles/nmat4998

First-principles molecular dynamics.

The behaviour of zeolitic imidazolate frameworks as a function of temperature was studied by means of density functional theory (DFT)-based molecular dynamics (MD) simulations, using the Quickstep module43 of the CP2K software package44. We used the hybrid Gaussian and plane-wave method GPW as implemented in CP2K. The simulations were performed in the constant-volume (N, V, T) ensemble with fixed size and shape of the unit cell. A timestep of 0.5 fs was used in the MD runs and the temperature was controlled by velocity rescaling45. MD simulations consisted of an equilibration period of 5 ps and production period between 60 ps and 220 ps, depending on the temperature, to study the dynamics of the liquid.

The exchange–correlation energy was evaluated in the Perdew–Burke–Ernzerhof (PBE) approximation46, and the dispersion interactions were treated at the DFT-D3 level47. The Quickstep module uses a multi-grid system to map the basis functions onto. We kept the default number of four different grids but chose a relatively high plane-wave cut-off for the electronic density to be 600 Ry, as already used in ref. 48, and a relative cut-off (keyword REL_CUTOFF in CP2K) of 40 Ry for high accuracy. Valence electrons were described by double-zeta valence polarized basis sets and norm-conserving Goedecker–Teter–Hutter49 pseudopotentials, all adapted for PBE (DZVP-GTH-PBE) for H, C and N or optimized for solids (DZVP-MOLOPT-SR-GTH) in the case of Zn. Representative input files for the molecular dynamics simulations are available as supporting information, and online in our data repository at https://github.com/fxcoudert/citable-data.

The unit cell studied for ZIF-4 (space group Pbca) is the orthorhombic primitive unit cell, which contains 272 atoms, with cell parameters a = 15.423 Å, b = 15.404 Å, c = 18.438 Å, and α = β = γ = 90°. These unit cell parameters were the result of a full geometry optimization (atomic and unit cell parameters) of the ZIF-4 crystal at 0 K. The density obtained is 1.25 g cm−3, in good agreement with the crystallographic density of 1.22 g cm−3 at room temperature, and all simulations were performed with this density. To test the validity of the constant-volume approach (dictated by reasons of computation cost), we ran shorter constant-pressure (N, P, T) simulations at 300 K and 1,500 K. These simulations both yielded an average density of 1.22 g cm−3 (with instantaneous fluctuations of ±0.09 g cm−3), showing that the variations of density upon heating and melting are minimal. Finally, constant-volume FPMD simulations at 300 K and 1,500 K were performed at a density of 1.6 g cm−3, and their structure factors are compared to those obtained at 1.25 g cm−3 inSupplementary Fig. 16.

12. First-principles molecular dynamics investigation on Na3AlF6 molten salt

• https://doi.org/10.1016/j.jfluchem.2016.03.004

2. Computational methods

2.1. Details of first-principles molecular dynamics

In this work, we employed the combination of IPMD and FPMD to improve the calculation efficiency. Initial configuration of ions system for molecular dynamics was prepared by packing ions randomly into a given simulation cell using the Packmol code 25(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0125). We consider Na3AlF6 molten salt as the industrial mole component of NaF 75%, AlF3 25%, so these simulation cells were composed of 60 Na atoms, 20 Al atoms, and 120 F atoms (a total of 200 atoms). The models of molten salt were firstly equilibrated with IPMD and then FPMD calculation initiated from the final converged structure of IPMD calculations. The method of starting FPMD calculation from the liquid prepared by IPMD has been employed successfully in the literature 26(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0130). Note that FPMD is quite insensitive to the quality of IPMD final converged structure, so the accuracy of potential parameters used in IPMD simulation don’t affect the final result of FPMD calculation. The IMPD simulations were run with the LAMMPS 27(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0135) software using the Buckingham potentials for Na3AlF6 molten salt system, and the potential parameters were referred from the literature 3(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0015). Verlet Leap-Frog algorithm was used with a time step of 1 fs to solve the equation of Newton motions. Ewald sums were used for all coulomb and multipolar interactions with a buffer width of 0.5 Å and an accuracy of 10−5 kcal/mol 10(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0050), 11(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0055), 12(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0060). The short-range interaction cutoff was set to 15 Å. Formal charges were used for Na (+1), Al (+3), and F (−1). The periodic boundary conditions were also applied on all sides of the model boxes to create an infinite system with no boundaries, so that the calculated results would be more convincible. To mix the system completely and eliminate the effect of initial distributions, the box systems were heated up to 4000 K in an NPT ensemble for 100 ps at 1.01 MPa, which means these simulations were run keeping the number of particles (N), pressure (P) and temperature (T) of the system constant. Then the hot liquids were cooled down at a rate of 1 K/ps to the melting point 1283 K. Another equilibrium at NPT ensemble for 100 ps was performed for relaxation. After these runs, convergence was achieved and the density difference between the initial and the final state fell below 1%.The resulting structure and velocities were used to start the following FPMD simulation.

FPMD simulation for Na3AlF6 molten salt has been launched with the CASTEP (Cambridge serial total energy package) code 28(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0140), 29(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0145). Perdew–Burke–Ernzerhof (PBE) exchange-correlation function was implemented in the generalized gradient approximation (GGA) 30(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0150). Ultrasoft pseudo potentials (USPP) have been employed for all the ion-electron interactions. The ionic cores are represented by USPP for Na, Al and F atom. The Na 2s22p63s1 electrons, Al 3s23p1 electrons and F 2s22p5 electrons are explicitly regarded as valence electrons. Dispersion force was also included by using the semi-empirical DFT-D2 method 19(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0095), 26(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0130), 31(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0155). Energy cutoff of 350 eV and a 1 × 1 × 1 k-point mesh was chosen for FPMD simulation. Time step for FPMD simulation was chosen at 1 fs to insure an energy drift less than 1 meV/atom/ps. FPMD simulation was done in a statistical ensemble with fixed particle number, volume and temperature (NVT) using the Nosé-thermostat 32(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0160). The simulation temperature (1283 K, melting point of Na3AlF6) was similar to IPMD, and the density of simulation model was set to the experimental value 33(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0165), 34(https://www.sciencedirect.com/science/article/pii/S0022113916300562#bib0170) 2.095 g/cm3. Periodic boundary condition was also employed in FPMD and simulation represent an infinite molten salt system. After FPMD being launched, the average pressure and resulting pressure in our NVT is 0.00014 and −0.003 GPa, respectively, which means the box volume for NVT calculation is reasonable. Finally, the trajectories of all atoms were collected for the following statistical calculations of structure and transport properties.

13. Ab initio mechanism revealing for tricalcium silicate dissolution

• https://www.nature.com/articles/s41467-022-28932-2

Methods

Atomistic model

Because the Ca3SiO5 (111) surface has been studied to have a greater possibility to form in practice25,37, we took the pure Ca3SiO5 (111) surface from our previous work25, which was cleaved from the M3 polymorph of Ca3SiO5 (obtained from CCSD38), as the M3 polymorph is the most frequently observed in industrial clinkers together with the M1 form39. The details of the DFT-based geometry optimization of the bulk crystal and surface slab were listed in Supplementary Methods. Considering the adsorption of water molecules on the Ca3SiO5 surface occurs even before contacting bulk water40, we firstly adsorbed isolated water molecules to saturate the dangling bond on the Ca3SiO5 surface (Fig. 1b). Then, we put a 20 Å thick layer of water with density of 1 g/cm3 on the Ca3SiO5 (111) surface (totally 333 atoms) for the dissolution simulations. The lattice parameters were 14.21 Å × 11.72 Å × 36 Å after setting a vacuum of 15 Å along z direction. After the detachment of the Caβ from its initial position to the state of D, we further constructed a new model (14.21 Å × 11.72 Å × 48 Å) by adding another 10 Å thick layer of water on the previous system (totally 498 atoms) to calculate the structural and dynamic properties of the equilibrium state.

AIMD simulations

All the AIMD simulations reported in this work were performed within the framework of DFT with the generalized gradient approximation (GGA) using the Perdew-Burke-Ernzerhof (PBE)41 functional and Grimme D3 correction42, which was implemented in the CP2K/Quickstep code43. The Core electrons were described by Goedecker-Teter-Hutter (GTH) pseudopotentials44,45 and the valence electrons were described by a mixed Gaussian and plane waves basis (GPW)46. The wave functions were expanded on a double-ζ valence polarized (DZVP) basis set along with an auxiliary plane wave basis set at a cutoff energy of 500 Ry. The Brillouin zone was sampled by the gamma approximation. During AIMD, the nuclei were treated within the Born–Oppenheimer approximation with a timestep of 0.5 fs for equilibrium simulation, while 1 fs for metadynamics simulations with the replacement of hydrogen by deuterium to accelerate the structural evolution without energy drifts31,47. The temperature was maintained at 300 K using a Nosé-Hoover thermostat48,49 coupled to the system with a time constant of 1000 fs in the Canonical ensemble (NVT). The convergence criterion for energy was set to 10−12 Hartree and for self-consistent field was 10−6 Hartree. All the system were first optimized to a stable state and then thermalized for at least 2.5 ps before the production run for statistical analysis. The production times for different simulations are shown in Table 1.

14. Solid-to-liquid phase transitions of sub-nanometer clusters enhance chemical transformation

• https://www.nature.com/articles/s41467-019-13509-3

Methods

Computational models

The structures of bare Au13, and Au8 supported on clean MgO(001) were constructed to study the structural dynamics of the clusters and its effects on catalysis. The Au13 cluster was initially built with a highly symmetric cuboctahedral (Oh) structure. O2 prefers to adsorb on the hollow site (Supplementary Fig. 1a), and H2O prefers to adsorb on the top site (Supplementary Fig. 1b). The Au13 cluster was simulated in a cubic cell of 15 × 15 × 15 Å3. A four-layer MgO(001)-p(4 × 4) slab was used for the support of the Au8 cluster. The cells were modeled under 3D periodic boundary conditions, and the slabs and their images were separated by vacuum with a length of 15 Å. O2 molecule prefers to adsorb at the interface between the cluster and the support (Supplementary Fig. 1c).

Density functional theory (DFT) calculation

The AIMD simulations were carried out using the freely available program package CP2K/Quickstep42 The DFT implementation is based on a hybrid Gaussian plane wave (GPW) scheme, the orbitals are described by an atom centered Gaussian-type basis set, and an auxiliary plane wave basis set is used to re-expand the electron density in the reciprocal space. Perdew-Burke-Ernzerhof (PBE) functional43 with Grimme’s dispersion correction44 was used. The core electrons were represented by analytic Goedecker-Teter-Hutter (GTH) pseudopotentials45,46. For valence electrons (5d106s1 for Au, 2s22p4 for O, 2s22p63s2 for Mg, 1s1 for H), the Gaussian basis sets were double-ζ basis functions with one set of polarization functions (DZVP)47. We performed the spin-polarized DFT calculations on Au13 and Au8/MgO. The spin state of O2 on Au13 and Au8/MgO is doublet and singlet, respectively, and the spin state of H2O on Au13 is doublet.

Free energy calculation

In the work, we have calculated the reaction free energy profiles by combining AIMD, constrained MD and thermodynamic integration. The reaction free energies (barriers) can be obtained by integrating potentials of mean force (PMF) over a chosen reaction coordinate48, i.e., the O–O distance of O2 dissociation reaction and the H–O distance for H2O dissociation. The PMF is calculated using a Lagrange multiplier method, by averaging the force applied on the system to keep the reaction coordinate constant in AIMD runs49. In the AIMD simulations, canonical ensemble (NVT) conditions were imposed by a Nose-Hoover thermostat under various temperatures. The MD time step is set to 0.5 fs. In AIMD runs, the trajectories of the first 5–10 ps were regarded as equilibrium periods to ensure equilibria of the systems, followed by another 5–15 ps of production periods for data analysis. For the size of systems studied in this work, the time scale of ~10 ps is sufficient to obtain well converged PMFs. as evident by the time accumulating averages shown in Supplementary Fig. 3a, b. Integrating the forces against the distance gives the free energy profile of the O2 dissociation reaction (Supplementary Fig. 3c). In the force-distance curves, three points crossing the force zero correspond to the reactant, transition and product state, respectively. Note that the hysteresis in our thermodynamic integration is very small. As shown in Supplementary Fig. 3d, when calculating the mean force at the O–O distance of 2 Å, we have started the structure models in both forward and backward directions, i.e., the models with O–O distances of 1.9 Å and 2.1 Å, respectively. The corresponding averages of the mean forces are very similar, i.e., −0.113 eV Å−1 for the initial structure with O–O distance to be 1.9 Å (blue) and −0.106 eV Å−1 for O–O distance of 2.1 Å (red), indicating good convergence of our PMF calculations. To calculate the entropy change, we firstly fit the temperature dependent free energy curves, and then the entropy was obtained by taking the temperature derivative. Similarly, the heat capacity was obtained by taking the temperature derivative of the fitted canonical caloric curve 〈E〉(T).

For comparison, the Gaussian 09 code50 was also used for geometry optimization calculation of O2 dissociation on Au13. The PBE1PBE51 functional was employed to optimize the geometries of Au-O2 complex. The aug-cc-pvtz52 and Lanl2DZ53 basis sets were employed for O atom and Au atom, respectively. During the structure optimization, all the atoms (Au and O atoms) were allowed to relax. The vibrational frequency calculations were carried out to identify the stationary points and transition states (TS) with zero and one imaginary frequency, respectively. The intrinsic reaction coordinate (IRC) calculation was also performed to verify that the transition state connects correctly to the expected minima.

15. Two-dimensional monolayer salt nanostructures can spontaneously aggregate rather than dissolve in dilute aqueous solutions

• https://www.nature.com/articles/s41467-021-25938-0

Methods

Classical molecular dynamics simulations

The classical MD simulations were carried out using the GROMACS 4.5 package59. The modeled systems consisted of aqueous ion solutions confined between two smooth walls (nanoslits). The ions were modeled using the Charmm27, Amber03, and OPLSAA force fields, to ensure that the results of the simulations were independent of the model employed. Water molecules were represented by the TIP4P53, TIP3P56, and SPC/E57 models. The Lennard–Jones (L–J) and electrostatic parameters of the ion and water models were listed in Supplementary Table 1. The solution–wall interactions were described by the L–J 10–4 potential function, corresponding to the integral of the L–J 12–6 potential of the graphene walls. Thus, a nanoslit with a width of 0.8 nm could accommodate one layer of aqueous solution at ambient pressure. Rigid graphene walls were also used to describe the ion–wall and water–wall interactions. The L–J parameters for carbon atoms were σc = 0.34 nm and εc = 0.3598 kJ/mol. Similar results from the simulations of solution confined between two graphene walls showed that the spontaneous formation of NaCl or LiCl nanostructures within confined water was not affected by the solution–wall interactions (as shown in Supplementary Fig. 10).

All MD simulations were performed in the constant lateral pressure and temperature (NPxyT) ensemble, with periodic boundary conditions in the lateral directions (x and y). The temperature and pressure were controlled by the Nosé–Hoover thermostat60,61 and Parrinello–Rahman barostat62, respectively. A cutoff of 1 nm was used for the L–J interactions, and long-range electrostatic interactions were treated by the slab-adapted Ewald sum method63.

Free energy calculations

Potential of mean force (PMF) profiles for different ion pairs were calculated using the umbrella sampling algorithm64. The force constant adopted for the harmonic bias potential was 5000 kJ/(mol·nm), due to the strong interactions between cation and anion. The harmonic force was used to constrain the ions at distances of 0.15–1.0 nm, in increments of 0.01 or 0.1 nm to enhance the resolution and smoothness of the PMF. Each distance interval was sampled for 10 ns, and the data obtained from the last 5 ns were analyzed using the weighted histogram analysis method (WHAM)65.

Ab initio molecular dynamics simulations

AIMD simulations were performed using the Quickstep module implemented in the CP2K package58. Ion–valence electron interactions were represented by Goedecker–Teter–Hutter (GTH) pseudopotentials66,67. The GTH-valence double-zeta-polarized Gaussian basis combined with a plane-wave basis set (with an energy cutoff of 280 Ry) was selected for the AIMD simulations. The Gaussian and augmented plane wave (GAPW) scheme was applied to obtain well-converged forces for the Na+ ion47,68. The BLYP exchange–correlation functional was used together with the Grimme dispersion correction (D3)69. A time step of 0.5 fs was used to ensure sufficient energy conservation for highly confined water systems. The temperature was maintained at 300 K in constant-temperature and constant-volume (NVT) AIMD simulations.

16. Revealing the intrinsic nature of the mid-gap defects in amorphous Ge2Sb2Te5

• https://www.nature.com/articles/s41467-019-10980-w

Methods

Molecular dynamics simulations

Thirty model structures of amorphous 225GST were generated using classical MD with periodic boundary conditions. The total number of atoms in each periodic cell was 315, with the individual number of atoms being Ge = 70, Sb = 70 and Te = 175. Initially, the atoms were placed randomly in a cubic simulation box, with the cell size calculated from the experimental density (5.88 g/cm3)57. A machine-learned interatomic potential for amorphous 225GST24 was employed for the classical MD simulations. The force field was obtained within the GAP framework58 by fitting a database of configurations, which were evaluated with DFT.

The MD simulations were performed using the LAMMPS code59 and the glass structures were generated following a melt-and-quench approach. The canonical ensemble (constant number of particles, volume and temperature, or NVT) was applied to keep the density of the simulated amorphous model close to the experimental value. The Nosé–Hoover thermostat, with a relaxation constant of 40 fs, was chosen to control the temperature fluctuations. A time-step of 1.0 fs was used to integrate the equations of motion during the MD simulations. For each glass model, the initial configuration was heated up at 3000 K with a 30 ps MD run to ensure that the system was melted. The molten structure was subsequently cooled down to 1200 K and equilibrated using first the NVT and then the NVE (microcanonical) ensembles with 30 and 10 ps MD runs, respectively, in order to obtain a well-equilibrated liquid structure at this temperature. The system was then linearly cooled down to 300 K at a quenching rate of −15 K/ps. At 300 K, the glass structure was equilibrated for 30 ps, followed by a 10 ps MD run with the NVE ensemble, to collect structural data. Finally, the temperature of the system was brought down to around 0 K at a quenching rate of −15 K/ps. The computational scheme used in this work corresponds to a total simulation time of 190 ps for each model.

An amorphous 225GST model of 900 atoms was also generated by GAP-MD simulation, following the same melt-and-quench protocol used for the 315-atom glass structures. In Supplementary Note 10 is discussed how the molecular dynamics simulations performed in this work can ensure the independence of the generated amorphous models.

Electronic-structure calculations

The output glass structures at 0 K from the melt-and-quench simulations were used as input configurations, for the 30 315-atom and the 900-atom amorphous 225GST models, to further optimise their geometry, using DFT implemented in the CP2K code, and to calculate the electronic structures60. The CP2K code uses a Gaussian basis set with an auxiliary plane-wave basis set61. All atomic species were represented using a double-ζ valence-polarised Gaussian basis set62 in conjunction with the Goedecker-Teter-Hutter pseudopotential63. The plane-wave energy cutoff was set to 400 Ry. The atomic geometry of each glass model was optimised, and the electronic structure was calculated by using the non-local PBE0 functional, with a cutoff radius of 3 Å for the truncated Coulomb operator64. The inclusion of the Hartree–Fock exchange provides an accurate description of the band gap and the localised defect states in our 225GST glass models, as was demonstrated in our previous work65. The computational cost of hybrid-functional calculations can be reduced by using the auxiliary density-matrix method66 (see Supplementary Note 11 for details), as employed in previous modelling studies of amorphous materials67,68. The Broyden–Fletcher–Goldfarb–Shanno algorithm was applied in the geometry optimisations and the forces on atoms were minimised to within 0.023 eV/Å. Electron trapping was modelled by injecting an extra electron in the relaxed ground-state glass structure and minimising the energy with respect to the atomic coordinates. The same non-local functional, as well as the same set-up, were used for the geometry optimisations of the glass systems with the extra electron.

17. Aging mechanisms in amorphous phase-change materials

• https://www.nature.com/articles/ncomms8467

Methods

Computational details

All the amorphous models contained 216 atoms and were generated with Perdew–Burke–Ernzerhof functionals based on the GGA33 using QUICKSTEP53, a mixed Gaussian and plane-wave code implemented in the CP2K suite of programs32. The Kohn–Sham orbitals were expanded in a Gaussian basis set with triple/double-zeta plus polarization quality, whereas the electron density was expanded in plane waves with a cutoff of 300 Ry. Scalar-relativistic Goedecker–Teter–Hutter pseudopotentials54 were employed. The Brillouin zone was sampled at the Gamma point only.

The amorphous GeTe models produced under different chemical substitution conditions were further annealed with the vdW-DF2 functional39 using PWSCF, a plane-wave code included in the Quantum Espresso package55. In this case, LDA norm-conserving pseudopotentials (Perdew and Zunger56) were used with a 34 Ry cutoff energy for the kinetic energy. The molecular dynamics simulations were performed with a 3.84 fs time step and a Berendsen thermostat (relaxation time equal to 2.7 ps). The absorption spectra were computed using the turbo-Lanczos approach48 with a total of 4,000 recursion coefficients.

Amorphous models preparation is detailed in the Supplementary Method 1.

Photothermal deflection spectroscopy

PDS is a highly sensitive technique to measure optical absorption (see Supplementary Method 2)57. PDS measurements on post-annealed and amorphous deposited phase-change films (sheet thickness d=200 nm) have been performed in cooperation with Reinhard Carius and Josef Klomfaß at the research center Jülich using a solid-state laser with a wavelength of λ=632 nm.

Impedance spectroscopy

Samples were prepared as a sandwich capacitor of 2 × 2 mm2 with multilayers (Au/Cr/Ti/TiN/GeTe/TiN/Ti/Au/Cr/glass substrate). The thickness of GeTe, the electrode (Au/Cr) and the antidiffusion layer (TiN/Ti) are 463, 105 and 27 nm, respectively. The electrodes were prepared by thermal evaporation and all other layers by DC magnetron sputter deposition.

The electrical impedance was measured via the Electrical Transport Option in a Quantum Design Physical Property Measurement System at the frequencies of 0.5–186 Hz from 4 to 150 K. Samples were annealed in the PPMS at different temperatures for one hour each with a ramp of 8 K.m−1. Results are analysed with the equivalent circuit of a single paralleled RC model.

18. Highly tunable β-relaxation enables the tailoring of crystallization in phase-change materials

• https://www.nature.com/articles/s41467-022-35005-x

AIMD simulations

Ab initio molecular dynamics (AIMD) simulations of Ge15Sb85 based on density functional theory (DFT) were carried out with the second generation Car-Parrinello-like MD scheme developed by Kuehne et al.68,69 and implemented in the CP2K package68,69. Two minimization steps were performed for each MD step. We employed the friction coefficients γL = 4.0 × 10−4 fs−1 and γD = 3.5 × 10−5 fs−1 (where γD corresponds to the intrinsic dissipation68,69, whereas γL is the coefficient of the Langevin thermostat). The Kohn–Sham wave functions were expanded in a triple-zeta plus polarization Gaussian-type basis set, and the charge density was expanded in plane waves with a cutoff 300 Ry. Scalar-relativistic Goedecker pseudopotentials and the standard exchange correlation potential parameterized by Perdew, Bruke, and Ernzerhof (PBE) were used70,71. The Brillouin zone was sampled at the Γ point. 1080 atoms were placed in a cell with dimension of 32.283 Å × 32.283 Å × 32.283 Å, and the corresponding number density was: 0.0321 Å−3 (ref. 72). The system was firstly randomized for 30 ps at 3000 K, followed by a subsequent equilibration for 30 ps at 1335 K. The melt-quenching model was obtained by quenching (3 K/ps) to 300 K. Note that this quenching rate is still rather faster than the experimental values. Re-annealed samples were obtained by equilibration at/re-heating to target temperatures (~458 K) for 500 k MD steps = 1 ns and 1 M MD steps = 2 ns. Finally, the system was re-quenched to ~285 K and equilibrated for 30 ps before further analysis. We computed the pair distribution functions g(r), and the structure factor S(q) was calculated from the partial distribution function by means of Fourier transform. The angular limited three-body correlation (ALTBC), which describes the probability of having a bond length r1 almost aligned with a bond length r2, was also computed from the AIMD trajectories (see Supplementary Fig. 13). More details can be found in ref. 72. We note that the increase in Peierls-like distortions is typically accompanied by increase in volume73. However, combining the efficient 2nd-generation Car-Parrinello scheme with constant-pressure simulations poses technical challenges. Thus, AIMD was carried out with a constant-volume protocol. The fact that we do observe a sharpening of the ALTBC peaks in spite of the constant-volume “constraint” indicates that the reinforcement of the Peierls-like distortions is a genuine effect.

19. The stability and electronic properties of novel three-dimensional graphene-MoS2 hybrid structure

• https://www.nature.com/articles/srep07007

Methods

In this work, two DFT packages have been used. The first-principles structure and energy calculations are mainly performed using the Vienna Ab Initio Simulation Package (VASP)43,44. Projector augmented-wave (PAW) pseudopotentials45 were used to account electron-ion interactions. The generalized gradient approximation (GGA) with the PBE function46 was used to treat the exchange-correlation interaction between electrons. The energy cutoff was set to 500 eV and 3 × 3 × 3 Monkhorst-Pack scheme was used to sample Brillouin zone. The full geometry optimizations are carried out with the convergence thresholds of 10−6 eV and 5 × 10−3 eV/Å for total energy and ionic force, respectively. DFT-D247 approach was used to describe the vdW interactions between layers.

The mixed Gaussian and plane-wave basis set code CP2K/Quickstep package48 has been used for the first principles molecular dynamics (FPMD) and Li diffusion barrier calculations. The wave functions of the valence electrons are expanded in terms of Gaussian functions with molecularly optimized double-zeta polarized basis sets (m-DZVP)49. For the auxiliary basis set of plane waves a 320 Ry cut-off is used. In the FPMD simulations the canonical ensemble was employed50,51 with a target temperature of 300 K, maintained with a Nosé-Hoover chain thermostat. Calculations of diffusion barrier of Li were performed using the climbing-image nudged-elastic-band method (CI-NEB)52 with nine points along the reaction paths.