AtomREM: Non-empirical seeker of the minimum energy escape paths on many-dimensional potential landscapes without coarse graining

TL;DR

AtomREM is an open-source, MPI-parallelized software tool that non-empirically finds minimum energy escape paths in high-dimensional potential landscapes, integrating with LAMMPS for diverse atomic systems.

Contribution

It introduces AtomREM, a novel parallelized solver implementing a non-empirical method for locating minimum energy paths without coarse graining.

Findings

Successfully applied to molecular systems

Effectively finds minimum energy escape paths

Open-source and compatible with LAMMPS

Abstract

Recently a non-empirical stochastic walker algorithm has been developed to search for the minimum-energy escape paths (MEP) from the minima of the potential surface [J. Phys. Soc. Jpn. 87, 063801 (2018); Physica A, 528, 121481 (2019)]. This method is based on the Master equation for the distribution function of the atomic configuration which has a nature to seek the MEP up along the valley of the potential surface. This paper introduces AtomREM (Atomistic Rare Event Manager), which is an MPI parallelized solver program package for executing this method, which yields minimum energy reaction pathways in terms of the microscopic evolution of atomic positions. It is open-source and released under the GNU General Public License (GPL). AtomREM interfaces with the LAMMPS Molecular Dynamics Simulator as a library of versatile potential functions for application to various systems. Examples of the applications to molecular and solid systems are presented.

Tables

Table 1. Table 1: External softwares

mt19937.f90	Subroutine generating random numbers by the Mersenne Twister method [20], distributed in http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/FORTRAN/fortran.html
libfwrapper.c	Wrapper for calling lammps from fortran code. This code is a modification of the one provided in LAMMPS (examples/COUPLE/fortran/)
USER-LAPLACIAN/	Package for calculating Laplacian to be installed in LAMMPS. This package is a modification of LAMMPS_LOCAL_HESSIAN package by S. Kadkhodaei [21], distributed in https://cmrl.lab.uic.edu/lammps_hessian.html.

Table 2. Table 2: Content of parameter file

Example	Explanation
&input	Start line for input file
Nstep = 12000,	Total number of time steps for simulation
Nwalker = 1600,	Number of walkers
Natoms = 7,	Number of atoms
temp = 0.1,	Temperature in the energy units of LAMMPS or analytical potential
tempFin = 0.4,	Final temperature in the energy units of LAMMPS or analytical potential
steptowrite = 50,	Number of steps to write the atomic coordinates and their energies
dt = 0.002,	Value of time step
ratio = 0.48	Value of initial delta function
mode = ”Initialization”	Name of mode = ”Initialization” or ”Langevin” or ”Reaction”
	Next data define the simulation box:
x_orig = -8.897999	x_orig, y_orig, y_orig are the coordinates
y_orig = -8.3346005	of the down left angle of the box
z_orig = -6.0167999
a_vec1 = 18.5668984 0.0 0.0	a_vec1, a_vec2, a_vec3 are the vectors,
a_vec2 = 0.0 17.8101998 0.0	which defines the supercell box
a_vec3 = 0.0 0.0 11.6429997
bounds = f f f	Boundaries of simulation box in the term of LAMMPS: f - non periodic, p - periodic
/	Finish line for input file

Table 3. Table 3: Format of output file ”path.dat”

Time step $\tau$	Average energy of atoms $E_q$ (from $q$)	Average energy of atoms $E_p$ (from $p$)	Force calculated for the average position of the walkers $F_{\mathrm{ave}}$	Number of walkers $N_w$
1	-0.39392953E+0001	-0.39392953E+0001	0.37627504E+0000	1600
50	-0.39394117E+0001	-0.39394117E+0001	0.34776722E+0000	1600
100	-0.39393155E+0001	-0.39393155E+0001	0.44461967E+0000	1600
150	-0.39389018E+0001	-0.39389018E+0001	0.64310856E+0000	1600
200	-0.39386394E+0001	-0.39386394E+0001	0.73788821E+0000	1600

Table 4. Table 4: Output files for mode ”Initialization”

File name	Columns in the file	Description
After_heating_energy.dat	$i_w{E}_i$333$i_w$ is number of walker444$E_i$ is an average energy of walker	Average atomic energy for each walker
Number_struct.dat	$i_w$   $U{n}_i$555$U{n}_i$ is a unique number of group   ${\tau }_{rlx}$666${\tau }_{rlx}$ is a time of relaxation	Information about groups and relaxation time
Structures.lammpstrj	-	Configuration data in the format of LAMMPS trajectory
Energies_of_structures.dat	$i_w$ $E_i$ $E_i^{\left(1\right)}$ $E_i^{\left(2\right)}$ … $E_i^{\left(N\right)}$ 777$E_i^{\left(1\right)}$, $E_i^{\left(2\right)}$, … , $E_i^{\left(N\right)}$ are the atomic energies for 1-st, 2-nd, …, N-th atom	Atomic energies and debagging information about all selected walkers
Chosen_walk_Structures .lammpstrj	-	Atomic configurations for chosen walkers
Chosen_walk_Number	$i_w$   $U{n}_i$   ${\tau }_{rlx}$	Information about groups and relaxation time for chosen walkers
Chosen_walk_Energy.dat	$i_w$  $E_i$  $E_i^{\left(1\right)}$  $E_i^{\left(2\right)}$ … $E_i^{\left(N\right)}$	Energies of each atom for chosen walkers
Chosen_BEST_walk.dat	-	All information about the slowest walker: number of walker and energies of each atom, number of group, time of relaxation, and coordinates of each atom. Last lines are utilizatble as atoms.dat for mode “Reaction”.

Table 5. Table 5: Output files for mode ”Reaction”

File name	Columns in the file	Description
atomic_coord_q.lammpstrj	-	Configuration data in the format of LAMMPS trajectory under biasing potential $V\left(𝒙\right)$
path.dat888please, see Table 3	$\tau$  $E_q$  $E_p$  $F_{min}$	Time evolution of the energy and the minimum force of the walker ensemble under biasing potential $V\left(𝒙\right)$
E_atoms.dat	$\tau$ $E_q$ $E_{\tau }^1$ $E_{\tau }^2$…$E_{\tau }^N$	Energies of each atom of the average position at time step $\tau$ under biasing potential $V\left(𝒙\right)$
saddle.dat	$typ{e}_i$ $x_i$ $y_i$ $z_i$999$typ{e}_i$ is a type of $i$-th atom, the $x_i$, $y_i$, $z_i$ are coordinates of $i$-th atom.	Atomic configuration of a saddle point before the relaxation stage
E_atoms_100XXX.dat	$\tau$ $E_q$ $E_{\tau }^1$ $E_{\tau }^2$…$E_{\tau }^N$	Energies of each atom for the XXX thread (walker) at time step $\tau$ at the relaxation stage
coordinates_rank_100XXX. lammpstrj	-	Atomic configurations for the XXX thread (walker) at the relaxation stage

Equations

$$d{x}_i=-\frac{{\partial }_{i}U(\mathbf{x})}{\Gamma }dt+\sqrt{\frac{2k_{\mathrm{B}}Tdt}{\Gamma }}{W}_i,$$

$${\partial }_{t}p(\mathbf{x},t)$$

$$\equiv$$

$$q(\mathbf{x},t)\equiv p(\mathbf{x},t)/\exp [-V(\mathbf{x})/k_{\mathrm{B}}T]$$

$$q(\mathbf{x},t+\tau )=\exp \{[{\hat{L}}_{\mathrm{Sm}}^{\prime }+{\hat{L}}_{\mathrm{rate}}]\tau \}q(\mathbf{x},t)$$

$${\hat{L}}_{\mathrm{Sm}}^{\prime }$$

$${\hat{L}}_{\mathrm{rate}}\equiv L_{\mathrm{rate}}(\mathbf{x},t)$$

$${\partial }_t\ln C(t)$$

$$F(\mathbf{x})$$

$$\exp \{[{\hat{L}}_{\mathrm{Sm}}^{\prime }+{\hat{L}}_{\mathrm{rate}}]\tau \}$$

$$\simeq \exp \{{\hat{L}}_{\mathrm{rate}}^{\prime }\tau /2\}\exp \{{\hat{L}}_{\mathrm{Sm}}\tau \}\exp \{{\hat{L}}_{\mathrm{rate}}^{\prime }\tau /2\}+O({\tau }^3).$$

$$\mathbf{x}(t+\tau )-\mathbf{x}(t)$$

$$=-\frac{1}{\Gamma }\nabla (U(\mathbf{x}(t))-2V(\mathbf{x}(t)))\tau +\sqrt{\frac{2k_{\mathrm{B}}T\tau }{\Gamma }}\mathbf{W}.$$

$$U(\mathbf{x})=\sum _{a<b}U({\mathbf{x}}^{(a)},{\mathbf{x}}^{(b)})+\sum _{a<b<c}U({\mathbf{x}}^{(a)},{\mathbf{x}}^{(b)},{\mathbf{x}}^{(c)})+\dots ,$$

$$E^{(a)}\equiv \sum _{b}U({\mathbf{x}}^{(a)},{\mathbf{x}}^{(b)})+\sum _{b<c}U({\mathbf{x}}^{(a)},{\mathbf{x}}^{(b)},{\mathbf{x}}^{(c)})+\dots (b,c,\dots \ne a)$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: AtomREM (Atomistic Rare-Event Manager)

Program Files: https://github.com/ryosuke-akashi/AtomREM

*Licensing provisions: * GNU General Public License 3 (GPL)

Programming language: Fortran 90 and C

Supplementary material: 1. Program code of AtomREM; 2. Results of simulation for argon solid and $(CH_{2})_{4}$ molecule

Nature of problem: AtomREM has been developed to help to find a reaction path on the atomic level using weighted Langevin mechanics on the inverse potential landscape. The method is describe the low temperature transformation of complex systems without artificial forces or/and collective variables.

Solution method: Recently [1] we designed a non-empirical scheme to search for the minimum-energy escape paths from the minima of the potential surface to unknown saddle points nearby with one dimensional application. The method is based on the Master equation and its solver algorithm is constructed to move the walkers up the surface through the potential valleys [2]. The stochastic algorithm uses a birth/death stochastic processes for numerous of walker in combination with the Langevin equation. Each walker obey the statistics of average movement on the reaction path under biasing potential. Under this consideration the reaction path is derived as the average of the walker distribution from stable to saddle states.

1 Introduction

Searching optimal values of bounded functions defined for many dimensional space concerns diverse scientific problems. In the context of physics (and also chemistry and biology), minimization of the potential energy with respect to the atomic positions is the representative of such problems–the optimal atomic configurations are related to (meta)stable states of molecules and solid systems. A closely related important problem is to locate the paths connecting two (meta)stable states with minimum potential energy barrier, as it specifies the possible chemical reactions and, more generally, structural deformations, as well as its approximate rate of occurrence as represented by the Arrhenius equation [1] . A situation where the destinations of the paths are unknown is especially interesting since it concerns seeking of reactions and deformations to form unknown products, which is the main focus of the present article.

Most existing methods to seek such reaction pathways using molecular dynamics and mechanics use collective variables and/or artificial forces; namely, enhance the movement of the atoms in certain targeted directions so that the desired reaction proceeds [2, 3, 4, 5, 6, 7, 8, 9]. However, those approach requires prior knowledge or assumption of the products to designate the collective variables or artificial forces and the solutions are inevitably subject to the human bias.

Apart from this mainstream, we have proposed a method based on the Langevin mechanics to seek the paths from known to unknown minima with minimum potential barriers without the collective variables or artificial forces [10]. This method has been demonstrated to be useful for seeking the reaction paths of the argon clusters with several tens of particles [11], which shows its potential broad applicability to systematic search for unknown reaction pathways described by atomic configurations. The purpose of the present article is to share the current version of the code package for this method as “AtomREM [12]”–Atomistic Rare Event Manager. The current version provides the following functionalities:

Nonempirical simulation of the paths from the known to unknown local minima with practical computational time up to at least several tens of atoms

• 2.

Output of the paths visualizable with open-source softwares such as OVITO [13, 14].

• 3.

Flat-MPI parallel execution for acceleration of the calculation

• 4.

Use of LAMMPS [15] as an external library for various potential functions.

In Sec. 2 we briefly explain the background theory and methods implemented in the current version of AtomREM published previously in Refs. [10, 11]. The contents and usage of AtomREM is described in Secs. 3 and 4, respectively. Some demonstrations are provided in Sec. 5. Section 6 is devoted to concluding remarks. Summary of the source codes and detailed discussion of the numerical behavior are also provided in Appendices for future developers.

2 Theory and method

Suppose any potential surface $U({\bm{x}})$ is defined on the space of configuration of $N$ atoms [${\bm{x}}\in{\bf R}^{3N}$; ${\bm{x}}=(x_{1},x_{2},\dots,x_{3N})$]. The task of AtomREM is to generate continuous trajectories non-empirically from a known local minimum of $U({\bm{x}})$ to unknown minima across the saddle points, on the basis of the overdamped Langevin equation

[TABLE]

where $x_{i}$ and $i$ are the coordinate and the index for the degrees of freedom, ${\bm{W}}=(W_{1},W_{2},\dots,W_{3N})$ is the vector whose components are randomly generated from the standard normal distribution at each step, $\Gamma$ is the friction constant, $dt$ is a timestep. At low temperature, the dominant rare trajectories across the saddle points nearby hardly occurs with straightforward simulation of Eq. (1) because of the exponential dependence of the rate of occurrence $\propto{\rm exp}[-\delta E/T]$ with $\delta E$ being the potential barrier height. On the other hand, under high temperature the atomic configuration is randomly perturbed strongly and information of the trajectory is lost or even system could be melt. AtomREM manages the generation of monotonic trajectories across the saddle points utilizing a recently proposed stochastic process [10, 11] which has a nature of tracking up the minimum energy paths (MEP) to the saddle points.

Here we briefly explain the background theory on the method. The rarity of the events across the saddle points with Eq. (1) can also be interpreted with the equivalent Fokker-Planck (Smoluchowski) equation [16]

[TABLE]

This is the equation that describes the probability distribution $p({\bm{x}},t)$ of the variable ${\bm{x}}$ at time $t$ evolved by the Langevin equation Eq. (1). In this representation, the rarity of the saddle crossing is interpreted as small amplitude of $p({\bm{x}},t)$ at the saddle points. The rarity in the long-time simulation is inferred by the Boltzmann distribution $p_{\rm B}({\bm{x}})\propto{\rm exp}[-U({\bm{x}})/k_{\rm B}T]$ as the stationary solution of Eq, (3). One can therefore conceive that, with a localized initial condition $p({\bm{x}},t_{0})=\delta({\bm{x}}-{\bm{x}}_{0})$, the quantity $p({\bm{x}},t)/p_{\rm B}({\bm{x}})$ must have sizable components that spread toward high-energy regime.

On the basis of the above consideration, the authors have examined the behavior of the distribution function

[TABLE]

with $V({\bm{x}})\equiv(1-\eta)U({\bm{x}})$ and $\eta$ being the control parameter. In the one dimensional harmonic potential case ($U({\bm{x}})=U(x)\propto x^{2}$) [10], the following properties have been found with the initial condition $p(x,t_{0})=\delta(x-x_{0})$: (i) it has the Gaussian form and its width gradually spreads with t; (ii) if $0<\eta<1/2$, its center goes up the potential surface in a short time and comes back to the potential minimum, whose initial velocity is proportional to $x_{0}$. The authors have demonstrated that these properties can be utilized to track up the “valley line (Fig. 1)” of general potential surfaces from the potential minima [10, 11]. Namely, in many dimensions, if we set the initial condition $p({\bm{x}},t)=\delta({\bm{x}}-{\bm{x}}_{0})$ with ${\bm{x}}_{0}$ in the middle of the valley line, the distribution $q({\bm{x}},t)$ goes up along the line keeping the well-localized Gaussian form. One can therefore draw the valley line extending toward the saddle point by tracking the short time evolution of the many-dimensional distribution $q({\bm{x}},t)$. The authors have indeed shown the successful visualization of the trajectories connecting the potential minima through the saddle points for argon clusters [11].

In summary, AtomREM solves the Master equation of $q({\bm{x}},t)$

[TABLE]

with

[TABLE]

Here $C(t)$ is the normalization coefficient for keeping the norm of $q$ unity. The average of the function $f({\bm{x}})$ is defined as $\langle f\rangle_{q(t)}=\int d{\bm{x}}q({\bm{x}},t)f({\bm{x}})$. The dynamics of $q({\bm{x}},t)$ is recast to the time evolution of walkers, which respectively have ${\bm{x}}$ with different values and their assembly forms $q({\bm{x}},t)$. By the Suzuki-Trotter decomposition [17, 18] Eq. (5) can be implemented as the repetition of usual Langevin evolution on a modified potential surface and stochastic copying and removal of walkers; the following formula is adopted in AtomREM

[TABLE]

The operator ${\rm exp}\{\hat{L}_{\rm Sm}\tau\}$ corresponds to the independent time evolution of the walkers by

[TABLE]

Note that, compared with the original Langevin equation Eq. (1), the potential surface is modified from $U({\bm{x}})$ to $U({\bm{x}})-2V({\bm{x}})$ and, with $0<\eta<1/2$, the modified potential is proportional to the inverse of $U({\bm{x}})$. This is the main driving force that makes the walkers go up the potential surface. The operator ${\rm exp}\{\hat{L}^{\prime}_{\rm rate}\tau/2\}$ is a multiplication of a scalar function to $q({\bm{x}},t)$ and implemented as the copying and removal of the respective walkers with the probability depending on the magnitude of the scalar function. For this task, an algorithm keeping the number of walkers strictly constant [19] has been implemented in AtomREM [11]. Parallel execution of the whole evolution step [Eq. (10)] is available.

We have explained above how the dynamics of $q({\bm{x}},t)$ is executed with a given initial position ${\bm{x}}_{0}$. Another important task of AtomREM is to generate the initial points that connect to various saddle points with low energy barriers. For this purpose “heat and relax” method [11] is used. Generally speaking, the relaxation of ${\bm{x}}$ by the steepest descent method behaves like the two-step process; fast relaxation onto the valley line and slow relaxation along it. One can therefore expect that, among randomly distributed (heated) walkers, the walkers which come late across any small energy threshold are expected to be on the valley lines. AtomREM also provides the function to get the initial positions which are likely to connect to low-energy saddle points on this principle. A simple description of the procedure will be appended later, though more detailed information is available in Ref. [11].

3 Contents of the package

AtomREM provides three functions:

Langevin–solve the overdamped Langevin equation Eq. 1, mainly for the purpose of finding local minima of the given potential function $U({\bm{x}})$ 2. 2.

Initialization–find the initial positions to start the tracking of the MEP of $U({\bm{x}})$ 3. 3.

Reaction–track the MEP of $U({\bm{x}})$ from the given initial positions toward the saddle points

It is composed of two directories, where the source codes, interfaces, and sample inputs for those functions are contained.

lammps_pot/–package implemented for utilizing the potentials of LAMMPS package (recommended) 2. 2.

analytic_pot/–package implemented with in-house analytic formula for the potential (only the Lennerd-Jones potential with open boundary condition is available in the current version)

3.1 Dependencies on external softwares

For usage of the potential functions implemented in LAMMPS code, AtomREM requires the users to build LAMMPS as a static library (for details see Sec. 3.4 Basic build options in LAMMPS manual page [15]). Also, the external codes summarized in Table: 1 are modified and redistributed in AtomREM.

4 Usage

AtomREM has to be built and used in Linux system. This software has been verified to run in CentOS6 systems with Intel compiler. Although in other systems the regular running is not guaranteed as is, we believe that it can be made run regularly by small modification of Makefile.

4.1 Installation

Import the source code by git.

$ git clone https://github.com/ryosuke-akashi/AtomREM

Enter either of the directories, lammps_pot or analytic_pot, edit Makefile for the proper compilers and include and library paths. and enter the following.

$ make

If the compilation is properly finished the executable a.out is made.

For lammps_pot/, additional procedure is mandatory before building AtomREM to install the user package of calculating Laplacian for Eq. (9) in LAMMPS.

$cp -r USER-LAPLACIAN${lammps_source_directory}

$cd${lammps_source_directory}

$ make yes-user-laplacian

$ make mode=lib mpi

After those steps static library liblammps.a is obtained. The directory

${lammps_source_directory} must be specified in Makefile of AtomREM for linking this static library.

4.2 Input description

AtomREM executables requires three input files:

params.in–Input parameters for controlling the calculations written in the fortran namelist format (order insensitive). Imported as a standard input. The parameters are described in Table 2 2. 2.

atoms.dat–Initial atomic position ${\bm{x}}_{0}$ to generate the initial distribution $p({\bm{x}},t_{0})=\delta({\bm{x}}-{\bm{x}}_{0})$. The format is described in Fig. 2. 3. 3.

in.pair–(Mandatory only for lammps_pot/) Part of LAMMPS input file where the potential function is defined. The format is described in Fig. 3.

The LAMMPS potential file declared in in.pair must also be put in the working directory. For example, with the lower in.pair in Fig. 3, the potential file named “CH.airebo” is needed.

4.2.1 Parameters in params.in

Most important parameters of the simulation are presented in the file of parameters (Table 2), which includes the numbers of walkers and atoms, the time of the simulation by the number of steps and the duration of the time step, temperature and initial value of delta parameter, the names of the input and output files, and also the simulation box.

The specific parameters for each mode (”Initialization”, ”Langevin” and ”Reaction”) will be explained in the corresponding subsections in the section ”Using AtomREM”.

4.2.2 Internal parameters

There are some tuning parameters in the source code. Although those have been optimized empirically so that the users do not have to modify them, some of them will be referred to in Sec. 4.4 for custom usage.

4.3 Output description

The program has a lot of output files and the number and kind of these files depend on the simulation mode, but the format of files is similar and they have just three main types. The first one is the sequential list of the time step, the energy, the force, the number of walkers, etc., which is related to the time evolution. The example of output data for the reaction path is in Table 3. The ’E_atoms.dat’ file has the similar format: the 1st column is a time step; the 2nd column is the total energy divided by the number of atoms; the 3rd column is the energy of the 1st atom; the 4th column is of the 2nd, …etc. The second main type is the atomic coordinates in the format of LAMMPS trajectories (output files with extension *.lammpstrj), which can be processed by visualization tools such as OVITO [13, 14]. For example, the ’atomic_coord_q.lammpstrj’ file has the atomic coordinates for reaction pathway given by the averaging with respect to the walkers. More detailed description of such output files for each mode is presented later in Sec. 4.4. Finally, the code has some miscellaneous output files, such as intermediate ones and default output of LAMMPS. Some of them are also explained later.

4.4 Running AtomREM

The a.out file is executed by entering the following line

$ ./a.out < param.in > out

Specifying param.in as the standard input is mandatory, whereas the standard output file name is arbitrary. When the calculation regularly ends, the following line is written at the end of the standard output:

End of simulation.

Flat-MPI parallelization is also possible; indeed, the default Makefile is compatible with the intel MPI. A standard parallel execution using mpirun is, e.g., as follows:

$ mpirun -np 16 ./a.out < param.in > out

The situation where AtomREM is useful is that the users know a (an approximate) metastable structure realized with a given potential function $U({\bm{x}})$ and want to seek for unknown transition pathways beyond the potential barriers. To achieve this goal, the users can run AtomREM in the three modes: ”Langevin”, ”Initialization” and ”Reaction”. The ”Langevin” mode is used to locate an accurate local minimum from a given initial position. The ”Initialization” mode is used for finding the initial atomic configurations as the starting points to the reaction paths, which are used as the input of the later execution with the “Reaction” mode. The algorithm of the finding of the initial atomic configurations is described in work [11]. The ”Reaction” mode is used for seeking of the reaction path from the initial atomic configuration. In order to choose the mode, you have to write the name of a mode in the file of the parameters. Below, we append a detailed description of the three steps.

4.4.1 Mode ”Langevin”

The usual Langevin equation Eq. (1) is executed independently for each process. The minimum number of walkers equal the number of the processes. This mode is provided for finding any potential minima closely related to a given atomic configuration. By executing this with low temperature the system reaches to the minimum through the steepest descent path. Although the users of AtomREM is assumed to know at least one metastable configurations of the system before the use of the code, the precise location of the minimum, which is essential for the later modes, could depend on the potential function used. The purpose of this mode is to correct such dependence. The output files of this mode are only ’E_atoms_100XXX.dat’ and ’coordinates_rank_100XXX.lammpstrj’, where ’XXX’ is the index of the processes. This mode is also internally called in mode “Reaction” after reaching the saddle point, to let the walkers move to other metastable configurations beyond the saddle point.

4.4.2 Mode ”Initialization”

Using the “heat and relax” algorithm described in the paper [11], the mode ”Initialization” generates initial atomic configurations that are supposed to be well on the valley lines toward the saddle points, or, “entrances” to the reaction paths [Fig.4]: Entrance A is a state with the largest energy deviation to the $y$-direction, and entrance B is a state with the largest energy deviation to the $x$-direction. Note that one initial configuration in principle corresponds to one reaction path.

The task with this mode is to evolve the Nwalkers walkers with the Langevin equation Eq. (1) first under temperature temp for duration Nstep and later under temperature tempFin ($\ll$temp) to relax the walkers along their steepest descent paths. When the temperature is switched, the potential energies of the respective walkers $E_{i}(i=1,2,\dots,Nwalkers)$ are checked. Among the walkers, only those whose energies satisfy the criterion $E_{\rm cut}-\Delta E<E_{i}<E_{\rm cut}+\Delta E$ are kept and others are killed (Fig. 5). Afterwards, the time steps are measured for the remaining walkers to come down to an energy threshold $E_{\rm threshold}(<E_{\rm cut}-\Delta E)$, so that the slowly relaxing walkers are detected. The measured time steps and atomic configurations of the walkers are recorded when they reach to $E_{\rm threshold}$ for later data processing. The parameters for the criteria are set internally, according to the initial temperature temp in main.f90.111The default values are $E_{\rm cut}=E_{\rm ini}+1.5\times temp$, $\Delta E=0.25\times temp$, $E_{\rm threshold}=E_{\rm ini}+0.15\times temp$, with $E_{\rm ini}$ being the potential energy at the start of this mode. One is recommended to modify the definition of those parameters with observation of After_heating_energy.dat so that many walkers are within the range [$E_{\rm cut}-\Delta E,E_{\rm cut}+\Delta E$]

The mode ”Initialization” yields several output files to help the choice of initial atomic configurations. See Table 4 for details. The main information in the output files are (i) the indices of the walkers kept at the temperature switching, (ii) the numbers of time steps required for passing through $E_{\rm threshold}$, (iii) the atomic configurations at the step of reaching $E_{\rm threshold}$, and (iv) the indices of atoms whose per-atom energies (explained later) show the largest deviations from their initial values. The last information can be utilized to a clustering analysis for obtaining the entrance configurations to diverse reaction paths [11]. Generally, the classical potential function can be represented as the sum of two-body, three-body and more-body terms as

[TABLE]

where ${\bm{x}}^{(a)}$ is the position of the $a$th atom. This form allows one to define the per-atom energy of the $a$th atom by

[TABLE]

AtomREM records the atomic indices which show the three largest values of $E^{(a)}$ for each walker and attributes it to a unique number $Un_{i}$ (e.g., if the third, second, and fourth atoms show the largest deviation of $E^{(a)}$ in this order, $Un_{i}=30204$). Choosing the most slowly relaxing walkers for each group of walkers having common $Un_{i}$ can yield the entrance points to the different reaction paths.222Note that the proper grouping according to $Un_{i}$ generally requires consideration of the spatial symmetry of the system, though it works to some extent without such consideration. For flexibility we do not provide specific tools for such data analysis, though the atomic configuration of the most slowly relaxing walker of all is written out in file Chosen_BEST_walk.dat for convenience.

4.4.3 Mode ”Reaction”

In the mode ”Reaction” the walker distribution evolves according to the equation [Eq. (5)] under temperature temp in order to seek the reaction path to the saddle point from a given initial point. During the simulation, the Nwalker walkers are equally distributed to the respective processes. After the evolution of timestep Nstep is finished, all the walkers are gathered to their average position, the number of walkers per process is changed to unity, and afterwards the final relaxation with Eq. (1) is executed under temperature tempfin independently for each process. As a result, with Nstep being large enough for reaching near the saddle point, some walkers go beyond the saddle point and reach to other metastable points. We thus have the reaction path in the form of the time evolution of the walker average of the atomic positions. The related output files are described in Table 5. For example, plotting path.dat and several E_atoms_100XXX.dats at the same time yield the figure like Fig. 9 given later.

During the first part of the simulation, the “reset and pullback” operation [11] is automatically executed with an interval defined internally. Specifically, the continuous evolution with Eq. (5) yields gradual spreading of $q({\bm{x}},t)$ and the simulation eventually becomes subject to incorrect departure from the valley line. To mitigate this instability, the distribution is occasionally reset to the delta function and afterwards evolved with the biasing potential $V({\bm{x}})=0$ (= usual Langevin equation) in a short time to let the walkers back onto the valley line.

The parameter $\eta$ defining $V({\bm{x}})$ has a crucial effect on the stability of the simulation [10, 11]. The necessary condition for tracking up the potential slope is $\eta<1/2$. With too large $\eta$, the walkers quickly stop climbing up the potential slope and go down. Smaller $\eta$ drives the walkers faster and longer up the potential slope, whereas it tends to induce the numerical instability more frequently. To our empirical knowledge the best way is to start with $\eta$ as close to $1/2$ as possible (set by ratio in param.in) and reduce it if the walkers do not climb up the potential surface. The following internal algorithm is adopted: (i) Save the energy of the initial configuration $E_{ini}$ and start the simulation with $\eta$ given in param.in; (ii) make the simulation for $N_{\rm check}$ time steps101010By default we have set $N_{\rm check}=400$; (iii) check the walker average of the potential energy $E_{q}$; if $E_{q}>E_{ini}$ continue the simulation; else, decrease $\eta$ to amplify $V({\bm{x}})$ and start over the simulation again with the new $\eta$.

5 Examples of simulations and validation

We have published the examples of the application of a previous private version of AtomREM to the surface reactions of Ar7, Ar13, and Ar38 in the paper [11]. Here we will show other types of transitions in the cyclobutane molecule and argon crystal. The example includes incorrect destruction of the molecule, the C-C bond-breaking transformation and collective sliding transitions as a elementary process of the solid-solid transition. The final example especially represents the advantage of the current method that generates the collective transformation without a priori assumption to the reaction coordinate as it is very difficult for solids to choose the collective variable [23].

5.1 Cyclobutane

We generated the reaction paths of cyclobutane: to butylene-like chain (A in Fig. 6), butadiene (B in Fig. 6) and tetramethylene (not shown).

In Fig. 6 the first two snapshots for the A and B paths are the same, where the second is the structure of the saddle point. After reaching the saddle point, the simulation of each process was independently executed under very low temperature (tempfin=$10^{-4}$ eV) and yielded different paths through the steepest descent paths or Intrinsic Reaction Coordinate (IRC). The only 4 IRC were unique and have led to the initial cyclobutane state (line 1 in Fig. 7), to the butylene-like chain (A in Fig. 6), to butadiene (B in Fig. 6) and to tetramethylene (not shown).

The time evolution of the potential energies along these paths extracted from “path.dat” are shown in Fig. 7. During the reaction simulation, the lines are serrated due to “reset and pullback” steps (See Sec. 4.4.3). The reaction pathways are also in the Supplementary movies (see Supplementary materials).

If the initial positions are chosen incorrectly, the destruction of the molecule occurs during the simulation (Fig. 8) and energy jumps to the high value (line 5 in Fig. 7). Even if the initial position is correctly on the reaction path, incorrect departure from the path driven by the random fluctuation can occur when the temperature (temp) is huge. The peak at the 24000-time step for the lines 1-4 (Fig. 7) is thought to correspond to the latter case; thanks to the reset and pullback step the walkers come back to the proper reaction path.

Note that the target of the present demonstration is just to show the capability of the method to track the reaction path on a given potential surface, regardless of the accuracy of the potential surface itself. There are a lot of work devoted to the research of cyclobutane, for example [24, 25, 26, 27]. In work [24] the energy of dissociation was calculated by the composite method CBS-QB3 of quantum chemistry, and intrinsic reaction coordinate calculations have been performed at the B3LYP/6-31G(d) level to ensure that the computed transition states connect the desired reactants and products, whereas we have used the classical potential function AIREBO [22]. The energies of the metastable states consequently differ. Also, the positions of the two released hydrogen atoms in the butadiene case is thought to be an artifact of the AIREBO potential, though it could be associated to possible dissociation of the hydrogen molecule H2.

Fig. 8 demonstrates the incorrect reaction path with interesting features. The third snapshot shows some hydrogen atoms, which are too close to the carbon atoms because these H atoms drift up their local potential gradients. Due to the too short H-C distance, the energy of the molecule becomes very high (line 5 in Fig. 7). There are several possible reasons for this behavior: the first is incorrect initial atomic configuration, which does not lie on the reaction path, the second is the high temperature of simulation or high time step, which leads to the high fluctuations, the third is not enough number of walker, which makes the simulation subject to fluctuations of the average position of the walkers. Although the simulation finally settles to fragments of molecules (final snapshot) after the relaxation, the metastability of this specific atomic configuration is obviously the artifact of the model classical potential. This is a typical behavior of the simulation when it fails.

5.2 The HCP-FCC transition of argon crystal

The most interesting case is a simulation of solid-solid phase transition (PT) because this type of PT has a high energy barrier or/and it occurs under some special conditions [28]. We applied AtomREM to the PT in the argon crystal. The theory for argon crystal [29] gives the two most stable structure for solid argon - $fcc$ and $hcp$, the authors of recent paper [30] performed the density-functional studies of argon at high pressure ($20-100~{}GPa$) and also have found that properties of $hcp$ structure are very close to the $fcc$. Fig. 9 shows the energy evolution and Fig. 10 and Fig. 11 demonstrate the side and top views of the atomic configurations during the simulation. For simulation we used the LAMMPS Lennard-Jones potential with $\epsilon=0.238~{}kcal/mol=0.0103207~{}eV$ [31], $hcp$ cell with $2\times 2\times 3$ size with 48 atoms. To find the initial configuration, we used the Langevin mechanics simulation in the mode ”Initialization” under temp$=0.5~{}K$ for 4000 time steps ($\tau=0.0005$) and with 1600 walkers. To determine the initial configuration for the later “Reaction” process, we compared the relaxation times among a group of walkers which exhibited two largest energy deviations for certain two neighboring atoms and selected the slowest among them (actually in the mode “Reaction” the path was initiated by those atoms, as indicated by white color in the second snapshot of Fig. 10). The simulations of the reaction pathways were performed at the mode ”Reaction” under the temperature $T=10^{-3}K$, $\tau=5\times 10^{-4}s$, 1728 walkers. The shape of the supercell was fixed during the simulation. To analyze atomic configurations and recognize the $hcp$, $fcc$ and $bcc$ structures, we used the OVITO tool (Open Visualization Tool [13, 14]) and adaptive common neighbor analysis with a variable cutoff [32].

We have found the saddle point related to the original $hcp$ structure and the layered $fcc-hcp$ structures (Fig. 9). To make the walkers overcome the saddle point efficiently, tempfin was set relatively high, due to which the plot of Fig. 9 showed fluctuation at the time step where $V({\bm{x}})$ was switched off. The energy barrier was $\sim 0.10/48=0.002083~{}eV/atom$ or $0.0480351089~{}kcal/mol$ for $hcp-fcc$ layered transition. This is reasonably smaller than the experimentally observed heat of fusion [33] of $0.2845~{}kcal/mol$.

A notable thing is that we realized the elementary processes of the solid-solid PT with large energy barriers without melting of the system. In the previous paper investigating the solid-solid PT in iron [28] the PT was induced by molecular dynamics with anisotropic compression so that the barrier height is lowered. AtomREM does not require such modification of the potential surface for realizing the elementary processes and enables us more straightforward comparisons to the experiments in arbitrary compression and strain setups.

A closely related method of the solid-solid PT has recently been demonstrated with the metadynamics [8], for example, in papers Refs. [23, 34]. There, the nucleation problems have been investigated using an entropy-related order parameter and enthalpy as the collective variables. In such an approach the information of detailed atomic positions in the reaction paths are coarse-grained and there remains the possibility that the simulation results crucially depend on the choice of collective variables. Our method, without loss of the detailed atomic positions, could also be useful for seeking optimal collective variables.

6 Conclusions and future work

This article has introduced AtomREM, which has been developed in order to study the possible structural transformations to the unknown potential minima without introducing the empirical collective variables and artificial forces. The code package is released under the GNU General Public License 3 and available in a public software repository [12] including example cases. Various potential functions are available via the interface to LAMMPS molecular dynamics package [15]. Demonstrations with the cyclobutane molecule and argon solid systems shows the performance of AtomREM.

Currently, AtomREM cannot be combined with the density functional theory calculations for more accurate potential surface. To achieve the comparable accuracy, combination with the Deep potential molecular dynamics method [35] is under way.

Acknowledgment We thank Taichi Kosugi for advice on the coding. This research was supported by MEXT as Exploratory Challenge on Post-K computer (Frontiers of Basic Science: Challenging the Limits). This research used computational resources of the K computer provided by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (Project ID:hp160257, hp170244, hp180184, hp190176).

7 Refences