Intrinsic ferroelectric switching in two dimension α-In${}_2$Se${}_3$

1 Introduction

In this notebook, we share the training database, force field model, training metadata, essential input files for density functional theory (DFT) calculations and molecular dynamics (MD) simulations, data analysis scripts, and selected original MD trajectories, as detailed in our paper [1]. Additionally, this notebook provides an automated workflow to reproduce the main results discussed in [1]. The model file for indium selenide, together with the complete training database and testing data, can be found in our GitHub repository. To access the repository, you can simply execute the git clone command (further details are provided in Section 3).

The directory is organized as illustrated in the following diagram:

• train: Contains the training dataset and the input.json file which holds the training metadata. Refer to Section 2.
• model: Stores the force field file frozen_model.pb.
• DFT: Includes a representative INCAR file for DFT SCF calculations that were used to construct training database. Refer to Section 2.3.
• test:
  • BondOrderAnalysis: Holds files analysing the local atomic environments represented in the training database. Refer to Section 4.1.
  • AvalancheDynamics: Offers selected MD trajectories and scripts for 1D domain walls driven by out-of-plane electric field in 2D α-In2Se3 analysis. Refer to Section 4.2.

2 Database Construction

2.1 Training database

The force field of In${}_2$Se${}_3$ utilized in this work is a deep neural network-based model potential, referred to as deep potential (DP). Details regarding the construction of the training database, DFT calculations, and metadata of the DP model were documented in our previous work [2]. Specifically, we adopted the DP-GEN, a concurrent learning procedure, to construct the training database (see details in Section 2.1). The initial training database contains DFT energies and atomic forces for structures derived from random perturbations of ground-state structures of $\alpha$,${\beta }_3^{\prime }$,${\beta }_1^{\prime }$,${\beta }^{\prime }$ ,${\beta }_2^{\prime }$,$\beta$ phases of In${}_2$Se${}_3$. The final training database comprises 22600 2D In${}_2$Se${}_3$ configurations and 2163 bulk configurations. You can access the training database in In2Se3/train.

2.2 DP-GEN

We employ the Deep Potential Generator (DP-GEN) to construct the training database. DP-GEN is a concurrent learning procedure consisting of three stages: labeling, training, and exploration, which together form a closed loop. Starting with an initial training database that contains DFT energies and forces for a few configurations, four DP models with distinct random initializations of neural networks are trained. In the exploration phase, one of these models is employed for MD simulations to explore the configuration space. Predictions (energies and atomic forces) are generated using all four models for each new configuration sampled from MD. For configurations that are well represented by the current training database, these four models should display nearly identical predictive accuracy. However, for those not well-represented, we expect the four models to produce scattered predictions with significant deviations. The maximum standard deviation of predictions from the four models serves as a criterion for labeling: configurations from MD with significant model deviation are labeled. The energies and atomic forces of these labeled configurations, as computed using DFT, are subsequently integrated into the training database for the next training cycle. Here, the maximum atomic force standard deviation, denoted as ε, is used as the labeling criterion. We introduce two thresholds, ε_lo and ε_hi; only configurations for which ε_lo < ε < ε_hi are labeled for DFT calculations. We set ε_lo = 0.12 and ε_hi = 0.25. The introduction of ε_hi is to handle the exceptions due to highly distorted configurations resulting from low-quality DP models (especially in the first few cycles of DP-GEN) or unconverged DFT calculations. The iteration stops when all configurations sampled from MD simulations satisfy a predefined accuracy across all four models. A primary advantage of the DP-GEN approach is its streamlined and largely autonomous data generation, minimizing human intervention.

2.3 DFT calculations

The first-principles DFT calculations performed by using the Vienna Ab initio Simulation (VASP) package. The projected augmented wave method is employed, and the generalized gradient approximation of the Perdew-Burke-Ernzerhof (PBE) type is chosen as the exchange-correlation functional. The energy cutoff is set at 700 eV, and the k-spacing is set at 0.3 Å^-1 . A sample INCAR file for the self-consistent field (SCF) calculations can be found in the In2Se3/DFT directory.

2.4 Deep Potential

The DP model, based on a deep neural network with the number of learnable parameters on the order of 10${}^6$, offers a robust mathematical structure to represent highly nonlinear and complex interatomic interactions while bypassing the need to handcraft descriptors that represent local atomic environments. Specifically, the DP model features a symmetry-preserving embedding network that maps an atom's local environment to inputs for a fitting neural network which then outputs the atomic energy; the sum of atomic energies yields the total energy. The original references to the DP model can be found in [3] and [4].

In this study, we utilized the smooth version of the DP model and employed the DEEPMD-KIT package for training. The cutoff radius is set to 6 Å, with smoothing starting at 0.5 Å. The embedding network follows a ResNet-like architecture with dimensions (25, 50, 100). The fitting network consists of three layers, each containing 240 nodes. The loss function is defined as:

$$\begin{array}{c}L(p_{ϵ},p_f,p_{\xi })=p_{ϵ}\Delta {ϵ}^2+\frac{p_f}{3N}\sum _i{\left|\Delta {\mathbf{F}}_i\right|}^2+\frac{p_{\xi }}{9}{\left|\Delta \xi \right|}^2\end{array}$$

Here, $\Delta$ represents the difference between DP predictions and training data, $N$ is the number of atoms, $ϵ$ is the energy per atom, and ${\mathbf{F}}_i$ is the atomic force of atom $i$. The prefactors $p_{ϵ}$, $p_f$, and $p_{\xi }$ are adjustable parameters. We increased $p_{ϵ}$ from 0.02 to 1, while reducing $p_f$ from 1000 to 1. The input.json file for training is located in the In2Se3/train directory.

2.5 Fitting perfomance

Here we compare the energies and atomic forces predicted by DFT and DP for all the structures in the final training database.(a) energy per atom and (b-d) atomic forces for configurations in the final training database. The insets provide the distribution of the absolute error.

3 Getting Started

This section explains how to directly run calculations on Bohrium. The model file for hafnia, together with the complete training database and testing data, is available in our GitHub repository. To access the repository, you can simply execute following commands:

Then you should also enter your email, password, and project_id into job.json file.

Utilize replace.sh to replace all the job.json files within the directory

4 Results

4.1 Bond Order Analysis

The Bond Order Analysis obtained with the Deep Potential Molecular Dynamics can be found in the directory In2Se3/test/BondOrderAnalysis. The following workflow reproduces Fig.S10 that analysis local atomic environments using Steinhardt’s bond orientational order parameters $Q_4$ and $Q_6$.

First, we employ dpmd to generate $Q_4$ and $Q_6$ values for entire train configurations. Here, We only show example for one configuration.The calculation process may take 1-3 minutes.

After the simulation done. Then, you can extract the $Q_4$ and $Q_6$ values from the dump_ql.xyz file(which is output file from above simulation) by following below line.

The generated file ele_q46.txt saves the $Q_l$ values for each atoms. By iterating the whole database through the above steps, you can obtain the $Q_l$ values for all configurations in the dataset. The entire process may take over 20 hours. After collecting the data, running the following script will display the distribution of $Q_l$ values.

You can run following script to show the $Q_l$ distrbution and the mean absolute force error of entire database.

After the calculation is completed, you will obtain a figure similar to the figure presented below.

Similarly, by calculating the $Q_l$ values for the trajectory structures of domain wall movement in DPMD by applying electric field, you can obtain the distribution of ql values during the domain wall reversal process. By running the following code at a temperature of 380K and an electric field strength of 2V/nm (the relevant text file in the MD_Q4Q6.zip),you can observe the $Q_l$ value distribution obtained during the domain wall flipping process within the overall $Q_l$ value distribution of the entire dataset.

After the calculation is completed, you will obtain a figure similar to the figure presented below.

4.2 Avalanche Dynamics

An example for calculating the domain wall movement can be found in the directory In2Se3/test/AvalancheDynamics/dpmd_run_script. You can also run the following workflow in the directory In2Se3/test/AvalancheDynamics/dpmd_run_script to compute the domain wall position during the reversal at 300 K in the presence of an out-of-plane electric field of 2 V/nm.

First, submit a pre-equilibration task to obtain the configuration after relaxation.

After the simulation done.To apply an out-of-plane electric field to the configuration obtained after pre-equilibration, submit a job following the instructions below and retrieve the molecular dynamics trajectory-"ele_efield0.2.lammpstrj".

After the simulation done. You can obtain the output trajectory file named ele_efield0.2.lammpstrj(This file is too large, exceeding the maximum limit for uploading files on GitHub). Following the steps below,you can obtain data for plotting the domain wall position as a function of time.

Then, combining the generated output file "'position_value.out' and 'position_err.out'". The combined position data located inIn2Se3/test/AvalancheDynamics/plot_data/dw_position

After the calculations complete, you will obtain a figure similar to the figure presented below.

An example for analysising the domain wall movement interface can be found in the directory In2Se3/test/AvalancheDynamics/post-analysis/dw_interface. Using the following command by selecting the frame at 564 ps from the trajectoryele_efield0.2.lammpstrj obtained as shown above to obtain the file for domain wall interface dz_dump5640_1.xsf visualization.

After the calculations complete, you can visualize the file dz_dump5640_1.xsf using VESTA to obtain the visualization shown above.