Practical Introduction to the MD Simulations of Ionic Systems

Chapter 12 Practical Introduction to the MD Simulations of Ionic Systems 12.1 Examples of MD Simulation of Ionic System For anyone interested in performing molecular dynamics simulations, there are many possible approaches. Nowadays, many kinds of MD programs and those for analysis are available. Choice of the method of the calculations depends on the objective of the researcher and the application of MD depends on scientific background of the researcher. Due to large number of possible selections, it may be difficult to choose the most suitable programs or conditions for starting the MD simulations. In such situations, fundamental steps taken in the study of simpler system with a minimum setting will be useful as a start for learning how to perform MD simulations. Then one can select the suitable methods for one’s purpose. Therefore, in Sect. 12.1, we will show some examples and carry out the exercise of MD simulations of ionic systems using a short program. Although the program is not optimized one for your machine, it is essentially the same one used in our studies and can be used for other purposes. In this exercise, we planned that one can run MD simulations and treat the data on “windows machine” by oneself. Some programs and examples will be found in Electronic Supplementary Material (ESM) (http://extras.springer.com) You will be able to find materials by searching the book’s ISBN. Please use contents of ESM at your own risk. In Sect. 12.4, fundamental setting in MD programs is explained. How to prepare initial configurations of melt, crystals and glasses are also explained with examples. After that, introduction of several products of software of MD simulations and for visualization available freely or commercially are given. For advanced features to acquire, one may need to consult with experts of related systems, manuals of programs. In principle, it should be avoided to use addresses of WEB sites in publication. However, especially in the present section, addresses of some WEB sites are included for convenience of users, in spite of the fact we cannot guarantee the validity of the address. © Springer International Publishing Switzerland 2017 J. Habasaki et al., Dynamics of Glassy, Crystalline and Liquid Ionic Conductors, Topics in Applied Physics 132, DOI 10.1007/978-3-319-42391-3_12 533 534 12 Practical Introduction to the MD Simulations of Ionic Systems Example 1: MD simulations of lithium metasilicate glass on “Windows” In the first example, we use a program for the lithium silicate glass that works on “Windows” with a personal computer. STEP 1 (Preparation of Linux like environment on Windows) Unfortunately, the situations of the readers (for machines, OS etc.) can be different case by case. Here we assume the case of the “Windows 8.1 or 10” machine with 64 bit, while most of MD programs work on the “Unix” or “Linux” platform with suitable compilers. Therefore, at first, linux-like environment is introduced on windows by “Cygwin” [1]. Cygwin user can skip this step; although addition of some options may be required later. Cygwin is a large collection of GNU and Open Source tools which provide functionality similar to a Linux distribution on Windows. Please visit the WEB page (Ref. [1]) and follow the instruction there. Please include GCC and Gfortran compilers. Additional installations of editors (such as vi) may be useful for modification of the input files, although one can also use the text editor on windows. This environment can be used for many purposes. Some knowledge for Linux or UNIX may be required, although usages of some commands (appeared after the prompt $) are included below. Caution may be required for the difference of code (for example, it is found in carriage return (CR) or line feed (LF)), in using different systems. • If you are a user of the Linux or Unix and Intel compilers, please try to use the file named “mdintel.exe” in the ESM instead of “lsomd.exe”. • If your environment is not suitable for using exe files found in the ESM, please skip step 1–4. STEP 2 (Preparation of files on your home directory) Put all files within “Example 1” folder in ESM in your home directory of Cygwin (You can do it by the usual operation on windows). By the default installation, the home directory may be located at C:\cygwin64\home\“User Name”. For starting the program, one can prepare some shortcuts during the installation. For example, click “Cygwin64Terminal” created on the desktop for starting Cygwin. Input the following command and then depress the enter key to see file list. $ls –al File names with information will appear. STEP 3 Structure of Files lsomd.exe is a binary code of the MD Program LSOMD for M2SiO3 systems. The test data used is for M ¼ Li. Input files Initial configuration of the system of 3456 atoms (Li 1152, Si 576, O 1728) is given in lslarsinf800ini.data. Please copy the configuration data to “lslarinfini.data” for the input of the program by the following procedure. 12.1 Examples of MD Simulation of Ionic System 535 For copying the data from A to B files, put the command “cp A B” after the prompt, $, then depress a return key. In this trial, please input the following. $cp lslarinf800ini.data lslarinfini.data The data included in this file are initial coordinates of ions and atoms, xn, yn, zn in Å. That is, x1, y1, z1, x2, y2, z2, ----x3456, y3456, z3456 followed by dx1,dy1,dz1,dx2, dy2,dz2-----dx3456, dy3456, dz3456. Here dx1 is for a difference of x1(t) and x1(t þ dt) and similar definitions are used for dy, dz for n ¼ 1, 2,---3456. Several parameters are necessary to be set in “lslarset.data”. “lslarset.data”: Setting of total steps, temperature, pressure, delta t, etc. Initial setting in the file is as follows. ¼¼¼¼¼ lslarset.data¼¼¼¼¼ 5000 100 33.15 800. 1.00E 15 0. 6.941 28.0806 15.9994 0.87 2.40 1.38 1.0155 0.8688 2.0474 0.07321 0.03285 0.17566 10.87 23.18 70.37 ¼¼¼¼¼¼¼¼¼¼¼¼¼¼¼¼ For the first line, number of steps, N, interval of steps, mitv, length of the cell (in Å), target temperature (in K), and step time (in s) and parameter for setting conditions. In the last figure, ‘0.’ means NVE condition, while ‘1.’ means NVT condition (using Gaussian thermostat). The value ‘target temperature’ does nothing when the NVE condition is chosen. In the following lines, potential parameters (see Sect. 9.1) for lithium metasilicate are given. In the second line, mass of Li, Si and O are given. In the third line, effective charges of Li, Si and O are given, respectively. In line 4, (5 and 6), potential parameters in Gilbert-Ida type for these species, ai, (bi, and ci) are given in the order of Li, Si and O in this case. Please distinguish integer and real figures carefully. Setting of these parameters depends on the problem to be examined as well as the properties of the system. For example, one needs to modify the setting considering following situations. Usually longer total time steps will be necessary for slow dynamics in melts, glasses, and crystals. The highest frequency accessible depends on the mint of the output. Shorter step time is required for the system with faster dynamics. If the system includes light atom such as hydrogen, it should be less than 0.5–1 fs. Step 4 Run the program Please run the program “lsomd.exe” using the lslarset.data. Configuration file lslarinfini.data is also used by the program. Put the following command. 536 12 Practical Introduction to the MD Simulations of Ionic Systems $./lsomd.exe <lslarset.data Here “./” is for using a current directory for execution of the program. After pressing the “enter key”, the program will start. At first the content of “lslarset.data” is repeated for confirmation. Please wait for a while. Then the status of the run with interval of mitv steps will be appeared on screen. Count, Temperature for Li, Si, O and total (in K), Pressure (in GPa), Energies (Total, Coulombic, Repulsive and for cicj term (in J/mol). The following is an example. --------------------------------------------------------------------------------------------------COUNT T_Li T_Si T_O Tmean P Etot Ecoul Erepul Ecc 10 822.208 808.818 804.311 811.027 0.9078E-02 0.7689Eþ07 0.8087Eþ07 0.2975Eþ07 0.2576Eþ07 --------------------------------------------------------------------------------------------------The execution time depends on your machine. You can see temperature and energy fluctuate during the run. Total energy of the system is maintained if NVE condition is chosen. Pressure of the system tends to show large fluctuations including negative values. Note that the pressure is sensitive to the setting of the box size and also the potential parameters used. The value has a relatively large error because it is obtained from a difference of large comparable values. OUTPUT FILES When a run stops, the following files will be found. lslar2traj.data: Coordinate of atoms obtained for every mitv steps are accumulated in this file. The data have similar sequences as lslarinfini.data but is repeated for N steps with interval of mitv steps. The data for dx, dy, dz are not included here. lslar3thermo.data: This output is the same one, which is shown on the screen during the run. “lslarinflast.data” is the last configuration for the input of continuing run. The structure of the file is the same as the “lslarinfini.data”. Pair correlation function, g(r) of the system was stored in lslar-gr.data. Please examine the size of output files and “real time” required for the calculation of 5000 steps in your system. If your file space is enough and run time is acceptable, further longer runs can be done by your own choice. Before starting the second run, copy the lslarinflast.data to lslarinfini.data. $./cp lslarinflast.data lslarinfini.data Please try further the longer runs for better statistics. During the following runs, output files are overwritten by new data. Such “overwriting” seems to be used in several MD programs and caution may be required if you need long successive runs. When you need old data, rename the old files. Probably it is better to prepare the script file for renaming or modify programs to add the run number. It may be also useful to add the name of the system, temperature, date, etc. to the file name with some rules of order and abbreviations. 12.1 Examples of MD Simulation of Ionic System 537 Fig. 12.1 Example of the plot of g(r) using the “lslar-gr.data” obtained by the 100,000 steps run (Li2SiO3 system at 800 K) using “Excel” STEP 5 Visualization of data After finishing run, one can plot the g(r) for each combination of ion pairs using the file named “lslar-gr.data” (or “lslar-gr.csv” if you skipped step 1–4) by several plotting programs. For example, the plot using “Microsoft Excel” after some setting of the appearance can be found in Fig. 12.1. In the “lslar-gr.data”, the first column of the data is a distance r. The values of g(r) are found from the 2nd to 7th column for Li–Li, Si–Si, O–O, Li–Si, Li–O and Si–O pairs, respectively. Please examine the characteristics of the structures. An example of the output file “lso.pdb” will be found in the ESM. The pdb file can be read from many visualization programs. One can use it to examine the structural details such as angles, distances, etc. Please try to use VESTA (Visualization for Electronic and STructural Analysis) [2] and open the pdb file. Select “Edit” tab and “bonds” tab and set Si–O bonds (connecting the atoms within 2.1 Å is useful in almost all cases) in the case of VESTA. Then select the “ball and stick” model. You can also use other styles such as “space filling model”, “polyhedra (after setting bonds)”. The file “lso.vesta” will be found in the same location, where the Si-O bonds are connected and polyhedra of SiO4 units are shown as in Fig. 12.2. 538 12 Practical Introduction to the MD Simulations of Ionic Systems Fig. 12.2 Examples of the visualization of the Li2SiO3 system (3456 particles) (near the glass transition temperature) at 800 K. (a) SiO4 units are shown by polyhedra. Li ions are shown with atomic radius to emphasize the positions and existence of ionic paths. (b) The same structure as (a) but with ionic radius of Li ion. (c) The same structure in ball and stick model. (d) The same structure in space filling model. Green: Li ions, Red: Oxygen atoms and Blue: Si atoms (and SiO4 units). These figures are visualized using the file “lso.vesta” (or “lso. pdb”) in ESM and the program “VESTA” [2] Please change the size of atoms or ions by yourself. You can also Zoom, Rotate and Move it. You can also measure distances or angles in the structures. Please follow the manual of VESTA or other programs used for further details. You can also visualize the structure obtained by the run you have performed. Please run the program inftoxyz.exe to change the “lslarinflast.data” to “lso.xyz” by the following command. ./inftoxyz.exe <lslarinflast.data The file named “lso.xyz” was created in the same directory. Many programs for visualization can read “xyz” format. Please modify the appearance by your own choice. STEP 6 Further analysis of the system Change temperature or box size gradually. How does the structure change? You can also examine more details by your own programs. 12.2 12.2 Example 2: Analysis of the Lévy Flight and Lévy (Alpha Stable) Distribution 539 Example 2: Analysis of the Lévy Flight and Lévy (Alpha Stable) Distribution 12.2.1 Example of Analysis for a Mobile Ion In this subsection, an example for the ionic motion in the lithium “metasilicate” is shown, and after that examples included in the ESM will be explained. It shows a Lévy distribution as found in the functional form of the self-part of the van Hove functions for all ions. This result shown here is another proof of the existence of Lévy flight with α < 2. Original time series for ionic motion is plotted in Fig. 12.3, where the absolute value of the displacement for Li ion jri(t)j in lithium metasilicate at 500 K during 1 ns runs is used. Here t is a sampling time measured by an interval of dt. (In this example, dt is taken as 4 ps). 12.2.2 Log Return of Data For examining the functional form of vibrational motions and jumps, one can use “log return” defined by following relation, which is often used in the analysis of financial time series (for example, see Ref. [3]). For the time series, jri(t)j(¼jrtj), log-return, Δlnjrtj is given by, Δln r t ¼ ln r tþ1 ln r t ;   r tþdt r tþdt rt rt  ¼ ln 1 þ jr t j jr t j Fig. 12.3 Time series of absolute value of the displacement for Li ion jri(t)j in lithium metasilicate at 500 K during 1 ns runs plotted against t. Here t is a sapling time measured by an interval of dt (in this example, dt is taken as 4 ps), plotted by using Mathematica 540 12 Practical Introduction to the MD Simulations of Ionic Systems Fig. 12.4 Log-return of the time series of absolute value of the displacement for Li ion shown in Fig. 12.3 Therefore, the log-return is an estimation of the expansion rate of jrtj. The distribution of jrtj is affected by the sudden jumps of the data, and the change of the absolute value of the jrtj is affected by the distance from the initial position, while the log-return does not depend on the absolute value of the jrtj. Log return of the time series in Fig. 12.3 is shown in Fig. 12.4. 12.2.3 Comparison of the Distribution of Log-Return and the Fitted Curve In Fig. 12.5, histogram for the probability distribution for the log-return (denoted by Z here) (orange) is shown, where a fit by the stable distribution is represented by the blue curve (blue). Thus the motion of ion is represented by the Lévy flight dynamics with α ¼ 1.46. 12.2.4 Further Analysis of Cumulative Distribution and Comparison of Lévy (Alpha Stable with α < 2) and Gaussian Distribution (Alpha Stable with α ¼ 2) To check the quality of the fit, cumulative distribution function is useful. In the long time scale, the distribution is affected by the exponential truncation of the tail. One should note that the exponential truncation (see Appendix A.2) in the selfpart of the van Hove functions is found in the log-log scale (For example, see Fig. 9.8 of Chap. 9), where the contribution is less than ~1 %, but is emphasized in the log-log scale. 12.2 Example 2: Analysis of the Lévy Flight and Lévy (Alpha Stable) Distribution 541 Fig. 12.5 Probability distribution, P, of the log return (denoted Z here) of the time series of absolute value of the displacement for Li ion shown in Fig. 12.3. Histogram (orange) is for the data obtained by MD, while the blue curve is for the fitted one with the stable distribution. Parameters for (α, β, γ, δ) is found to be (1.46, 0.039, 0.073, 0.028) by the maximum likelihood method. Here the tail found in longer distance was omitted Fig. 12.6 Comparison of the pattern of the cumulative distribution of the log-return (blue dots) and that of the fitted one for stable distribution (red curve) 1.0 0.8 0.6 0.4 0.2 0.4 0.2 0.2 0.4 In Fig. 12.6, the function of the log-return of the time series (blue dots) and the fitted one for stable distribution (red curve) by the maximum likelihood method [4] obtained by using Mathematica are shown. In Fig. 12.7, comparison of the stable α distribution with α ¼ 1.46 (red) and corresponding distribution with α ¼ 2.0 (blue dot) is shown. The difference is observed in both the tail and the sharpness of the function. 12.2.5 Free CDF Files with Manipulation We prepared free Computable Document Format (CDF) files using Mathematica [5] for treatment of the data in the folder named “Example 2” in the ESM. It works on Windows. Please copy the contents of the folder in your Documents directory (by usual procedures on your machine). 542 12 Practical Introduction to the MD Simulations of Ionic Systems Fig. 12.7 Comparison of the pattern of the distributions. Red curve is the fitted one with (α, β, γ, δ) ¼ (1.46, 0.039, 0.073, 0.028) and blue dashed one is α ¼ 2 (other parameters are set to be the same as red one) 4 3 2 1 0.4 0.2 0.2 0.4 If you don’t have Mathematica, you can use a free CDF-Player which can be available at Ref. [6]. There is a free CDF file named “Levy -stable alpha- distribution of type 1. cdf”. It includes interactive manipulation of four parameters in the stable distribution (see Appendix A.2 for details of the distribution). If the file does not work well in your environment, please try “Levy -stable alpha- distribution of type 1_9.cdf” files. You can check how the distribution changes with parameters by yourself. There are several examples of CDF files (Analyses ---.cdf) for time series, jri(t)j data of some arbitrary chosen particle obtained by MD runs, which were fitted to the stable distribution. You can treat the file using the CDF player. Please change the size of figures to see the details for the motions Sample 1: Jump motion for a Li ion in lithium “metasilicate” glass at 700 K. Lévy flight dynamics is clearly observed. Sample 2: For an EMIM ion in the ionic liquid, EMIM-NO3 at 300 K. The motion is not a clear jump but the Lévy (alpha stable) distribution is clearly observed. Sample 3: For the arbitrary chosen Li ion at 800 K in lithium “disilicate” glass. The ion is localized for a long time (The time series during 12 ns run is shown here.). In this case, both fittings to normal and Lévy stable distributions are compared. The distribution is Gaussian like except for the contribution of the peak due to the jump found at the last of the time series. Note that the numbers of the sampling points are not the same in each example. The file named “Simulation using stable distribution.cdf” (and/or that with_9) is for the simulation of the displacement of ion reproduced by the random sampling of 1000 points from the stable distribution with a characteristic index α chosen by you. When α ¼ 2, the motion is the Brownian motion. You can find strong heterogeneous motions with intermittency with a smaller (<2) value of the index α. You can compare the actual data and the randomly sampled data. The free CDF file named “MD-Soft-Core Model.cdf” is for the MD simulations of soft-core model (exponent n ¼ 12). Although this is not an ionic system, readers 12.4 Fundamental Usage of MD Programs 543 can have some idea of what is done in the MD simulations. In this example, the periodic boundary condition is not used and particles are bounced back from the wall. Except for this point, the program do the MD simulation by the Verlet algorithm with your choice of the number of particles, temperatures and system size. Color of particle changes when the velocity of the particle changes. 12.3 Example 3: Examining Movies Two examples of movies are included in the “Example 3” directory in the ESM. Since their visibility depend on your machine environments, we prepared two files with different formats. We hope that you can see at least one of them. The file named “EMIM-NO3 100 K-1.wmv” is for EMIM-NO3 glass at 100 K shown by the stick model. This can be open by the “media player”. Try to click the file name on Windows. In these movies, periodic boundary condition is used as usual. If the center of mass position of the ion moves to the next image cell, the ion disappeared and the corresponding ion comes in at the same time. Thus the number of ions is kept constant. The file named “disilicate crystal.swf” is for lithium disilicate crystal in the metastable form at 1400 K. The format is for adobe flash player. Please try to open it from Internet Explorer. If it does not work, please check if the adobe flash player was installed in your system or not, by visiting the WEB site of adobe (Ref. [7]) and follow the instruction there. It can be seen that how the mobility of atoms depends on the direction of the crystal. For explanation of other files, see readme.txt or Readme---.docx in each folder in the ESM. 12.4 Fundamental Usage of MD Programs In followings, we explain a general (and may be minimal) example for input and output in MD programs. 12.4.1 INPUT of the MD Programs If you use a program suitable for your purpose, each MD program has its own input and output format and may have examples to show how to use it. 544 12 Practical Introduction to the MD Simulations of Ionic Systems Typical input files in the works using MD simulations of glasses contain 1. Setting of potential parameters and functions 2. Setting of ensemble (NTP, NVE, etc.) and some options such as the shape of the MD cell. 3. Initial configuration, velocity and system size. 4. Setting of runs for analysis (step time and total number of steps at each temperature). In the case of MD simulation of super-cooled liquids and/or glasses, cooling schedule must be planned. 5. Several settings for cutoff values (for repulsive term and/or Ewald summation). 12.4.2 Preparation of Initial Configurations for Crystals Crystal structures can be used for checking the quality of the potential parameters, for the preparation of initial configuration for liquids or glasses, as well as examination of crystals themselves. Here we briefly explain how initial configuration for MD simulations of crystal is prepared as exemplified by lithium or sodium metasilicate crystals. These crystals are orthorhombic and space group is Cmc21 with Z ¼ 4 [8]. Li2SiO3 a0 ¼ 9.36; b0 ¼ 5.395, c0 ¼ 4.675 (in Å) Na2SiO3 a0 ¼ 10.52; b0 ¼ 6.07, c0 ¼ 4.825 (in Å) Space group: Cmc21, where the meaning of each notation is as follows. C: base-centered lattice. m: mirror planes c: glide planes (mirror planes involving reflection and a translation parallel to the plane). (The plane can be a, b, c, n or d) 21: screw axis Treatment of crystal structure is based on the “International Tables for Crystallography” volume A [9]. General coordinate of the present systems are given by, 1. 2. 3. 4. 5. 6. 7. 8. x,y,z x, y,z þ 1/2 x, y, z þ 1/2 x, y, z x þ 1/2, y þ 1/2, z x þ 1/2, y þ 1/2, z þ 1/2 x þ 1/2, y þ 1/2, z þ 1/2 x þ 1/2, y þ 1/2, z 12.4 Fundamental Usage of MD Programs 545 Some software can generate the crystal structure by choosing the space group and input of coordinate in the unit cell. As an alternative way to obtain the basic cell for MD runs, the atomic coordinate, rα, within an asymmetric unit is expanded to ri in the whole crystal lattice by using symmetry operators. rα • Pm þ Tm ¼ ri ; ð12:1Þ where Pm and Tm are (3  3) matrix of linear part and (3  1) column matrix (translational vector), respectively for general transformation. The linear part (orientation or length, or both) of the basis vector a, b, c can be transformed into1 new vector a0 , b0 , c0 by the (3  3) matrix Pm ¼ 0 P11 P12 P13 @ P21 P22 P23 A, that is, (a0 ,b0 ,c0 ) ¼ (a,b,c)Pm. This is combined with translaP31 P32 P33 tional transformation. For example, obtaining the coordinates x, y, z þ 1/2, from atomic coordinate xα, yα, zα, following operation Pm 0 1 ð xα ; yα ; z α Þ @ 0 0 0 1 0 Tm 1 0 1 1 0 xi 0 0 0 A þ @ 0 A ¼ @ yi A zi 1=2 1 ð12:2Þ is applied for the (8b) positions in the atomic coordinates shown in the next Table 12.1 [8]. Both matrix and general forms also can be found in the WEB site [10] In Fig. 12.8, primitive cell for the lithium metasilicate crystal is shown. Usually, basic unit of MD simulations is chosen to have equal (or comparable) axis lengths by repeating the crystal units. In this example of lithium metasilicate crystal, a supercell formed by repetition of axis length as a0  2, b0  3, c0  3 is taken as the basic cell of MD simulations. The cell contained 2  3  3 of Z ¼ 4, namely 72 of Si atoms, 144 of Li ions, and 216 O atoms, and the total number of atoms is 432. The space group of the system is Cmc21. The supercell thus prepared is shown in Fig. 12.9. This system size can be used to examine some characteristics of crystal structures, but might be too small for examining dynamics of them because of long length scale repetitions of the structures affect the dynamics. Table 12.1 Atomic coordinate of lithium (and sodium) metasilicate crystal [8] Atom Li(Na) Si O(1) O(2) Position (8b) (4a) (8b) (4b) X 0.160(0.166) 0.000(0.000) 0.141(0.130) 0.000(0.000) Y 0.320(0.339) 0.164(0.166) 0.321(0.286) 0.100(0.077) Z 0.000(0.000) 0.537(0.564) 0.450(0.500) 0.860(0.895) 546 12 Practical Introduction to the MD Simulations of Ionic Systems Fig. 12.8 Unit of crystal in Li2SiO3 system [8]. Pale blue, blue and red spheres are for Li, Si and O, respectively Fig. 12.9 Example of basic cell of MD for Li2SiO3 crystal. Supercell consists of 2  3  3 of primitive cells. Colors are the same as in Fig. 12.8 The required system size depends on the property to be examined and it should be judged by the purpose of the work. Especially, oscillation of the system is affected by the symmetry of environments as well as that of the inner structure of the crystal unit, and therefore caution is required in how to prepare the supercell. 12.4.3 Data Base for Crystal Structures In the case of crystals, one can use the data from data bases such as AMCSD (American Mineralogist Crystal Structure Database) [11], COD (Crystallography Open Database) [12], WWW-MINCRYST (Crystallographic and Crystallochemical Database for Minerals and their Structural Analogous) [13]. Several data format (in xyz, cif and pdb) are used for initial configurations and changes of the format may be necessary to adjust the MD code. 12.4 Fundamental Usage of MD Programs 547 Configurations in the basic cell of crystal structure can be prepared by using the information of the space group as already mentioned. It can be done by using software, such as, “CrystalMaker®” [14], which can create a supercell of n1  n2  n3. Usually, it seems convenient to set the length of three axes to be comparable in 3D system as already mentioned; otherwise, different setting for calculations of cut off length, Coulombic force by directions will be required. Some other cautions are required in the preparation of the super cell. 1. Number of repeating cell may affect the properties of the crystal especially the dynamics. 2. Some atoms may be automatically added at the edge of the crystal by a software. It is usually unnecessary, when one uses PBC. One needs to eliminate them. 3. Some crystals contain several sites for one atom, with occupation number being less than 1. This should be carefully treated when it is used for an initial configuration of MD. 12.4.4 Initial Configurations for Melts and Glasses For melts and glasses, it can be started from the crystals or random configurations equilibrated at high temperatures. One may need to check if the final configurations after equilibration are affected by the initial configuration or not, especially when using crystal structures. For example, in the case of silica, it is known that the long range network structure remains for a long time in the structures. This trend is more remarkable in the strong system than in fragile systems [15]. Of course, several kinds of molecular structure editors can be used for preparing and/or modifying structures. The system in the glassy state is obtained by cooling or by compressing the melt. Some details of the glassy state depend on the cooling schedule. Example of a configuration of lithium silicate glass (after equilibration) is found in the ESM, and such a structure can be used for the initial structures of other melts. Required system size depends on the properties and status of the system considerably (see Sect. 8.5.1 for more details). Generally, in the melt and glass, the structure has shorter length scale than the crystal. Therefore in the former with periodic boundary condition, relatively smaller system can be used, as long as the effect of ghost particle (image of the particle in the basic cell) in the surrounding cells of MD simulations is negligibly small. In the lithium silicate glass, if one fixed the particle within a wall (for one of the 6 faces) of the basic cell, it affects mobility of other particles. The length scale thus determined was about ~8 Å at around 700 K for lithium metasilicate and the system size of L ~ 17 Å (longer than twice the length scale) for 432 particles seems to be the minimum size for examining of structure and dynamics of the present system in the liquid and glassy states. The larger size is required for stronger (rich in network formers) systems even for the liquid and glass [16], and it also means necessity of 548 12 Practical Introduction to the MD Simulations of Ionic Systems longer equilibration time of the system. Of course, practical system size depends on the available resources for calculations and “real” times necessary to do simulations. For the initial configuration of melts or glasses, random configurations at high temperature (for example, ~3000 K in silicate) are used frequently. If one start from the random configurations, it may be better to remove or change the position of particles being too close to other particles in the beginning of runs. The problem depends on the characteristics of potential parameters used. In the beginning of equilibration of the system, several pre-relaxation methods, which may be included in the software, can be used. Smaller time steps or modification of mass will be helpful to suppress unstable motion in the beginning. MD programs may have examples of initial configurations, tutorials and exercises, which are useful to begin simulations with/without some modifications, although check of the quality of parameters or methods used is the responsibility of the user. 12.5 Output of the MD Programs Most fundamental output of MD is time series of particle positions, that is, trajectories. From these values, many secondary outputs can be drawn. Further analyses by the user also can be done. In typical output files, fundamental information such as temperature, energies, pressure, pair correlation functions are included and they will be useful for checking the setting and/or for further treatment. 12.6 Software for MD Simulations Many kinds of MD programs have been developed including commercial ones and researchers can choose the ones suitable for their own purpose. It also depends on the available potential parameters, the functional forms and the ensembles to be used. Some of them, which can be used for the study of ionics, are introduced here, although for details the user should consult the developer. Amber [17] is a collective name of program. It is designed for the biomaterials but can be used for the ionic materials such as ionic liquids. Several force fields such as, general amber force field [18] can be used for the simulation of ionic liquid. Of course, other potential parameters tuned for individual systems are also useful. The “sander” code was for the classical molecular dynamics for CPU, while “pmemd” code was for GPU and is included in the version 11. After the Amber 12, “sander” has been added to the “pmemd” coding for both CPU 12.7 Software for Visualization 549 and GPU. Official support for the Intel Xeon Phi architecture also started in Amber 14. DL_POLY is a general purpose serial and parallel molecular dynamics simulation package (DL_POLY_4.08 is the latest version (March 2016)) [19, 20], which can be used for ionic systems and several exercises and demonstrations for them can be found. The Car-Parrinello Molecular Dynamics (CPMD) [21] is a code designed for the ab initio MD simulations, where the parallelized plane wave/pseudopotential implementation of Density Functional Theory is used. Large-scale Atomic/Molecular Massively Parallel Simulator, LAMMPS [22] is a molecular dynamics package and LAMMPSCUDA is a molecular dynamics package, which can use GPU and have been applied for ionically conducting glasses [23]. The OCTA [24] [OCTA is Greek “8”, which is a 90 rotation of the mathematical symbol of infinity.] is flexible and expandable programs for meso-scale simulation and it seems to be suitable for the multiscale simulations of polymer etc. MODYLAS [25] is a general-purpose, molecular dynamics simulation program suited for the simulation of very large physical, chemical, and biological systems. In the program, the Coulombic term can be calculated by the multipole method. Many of MD programs have been coded by Fortran (The name “Fortran” tends to be used for after Fortran 90, although “FORTRAN” had been used before.) and/or C (Cþþ). Therefore, some knowledge of them may be required. Traditionally, the former has been used for long time because of suitableness of numerical processing and there are many cumulative programs related to it. Its method for memory allocation has said to have merits for parallel computation. Some comparison of Fortran and C has been done by several authors. For example, see Ref. [26]. 12.7 Software for Visualization Nowadays, there are many software including commercially available ones and free ones for visualization of structures obtained by MD. One can select the suitable software for any particular purpose. Some of them are introduced here. Some tools can treat the output of MD codes, directly. For example, the program VESTA [2] can be used for the visualization of many kinds of files including CONFIG file used in DL_POLY code. VMD (Visual Molecular Dynamics) [27] can be used for visualization of instantaneous structures as well as making movies (for example, using HISTORY file obtained by DL_POLY and/or other trajectory files). 550 12 Practical Introduction to the MD Simulations of Ionic Systems References 1. http://www.wolfram.com/index.php. Accessed 12 Aug 2016 2. K. Momma, F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011). http://jp-minerals.org/vesta/en/. The address was confirmed to be valid on Oct. 2015 3. C. Brooks, Introductory Econometrics for Finance, 2nd edn. (Cambridge University Press, 2008) 4. http://www.wolfram.com/cdf-player/. Accessed 12 Aug 2016 5. http://www.adobe.com/software/flash/about/. Accessed 12 Aug 2016 6. “International Tables for Crystallography” volume A. WEB version is available at http://it. iucr.org/. Accessed 12 Aug 2016 7. http://www.adobe.com/software/flash/about/. The address was confirmed to be valid on 12th Aug. 2016 8. V.H. Seemann, Acta Cryst. 9, 251 (1956) 9. “International Tables for Crystallography” volume A. WEB version is available at http://it.iucr.org/. The address was confirmed to be valid on 12th Aug. 2016 10. http://www.cryst.ehu.es/cgi-bin/cryst/programs/nph-getgen?list¼new&what¼gen& gnum¼36. In this address, the number 36 means Cmc21. Please put the number of the space group in that page (These URLs are conformed to be valid on 26th Dec. 2014.) 11. http://rruff.geo.arizona.edu/AMS/amcsd.php. The address was confirmed to be valid on 19th Feb. 2016 12. http://www.crystallography.net/. The address was confirmed to be valid on 19th Feb. 2016 13. http://database.iem.ac.ru/mincryst/index.php. The address was confirmed to be valid on 19th Feb. 2016 14. http://www.hulinks.co.jp/software/c-maker/. The address was confirmed to be valid on 19th Feb. 2016. After the version 9.2.2, CrystalMaker® can treat STL files, which are supported by many 3D printers 15. J. Horbach, W. Kob, K. Binder, C.A. Angell, Phys. Rev. E 54, R5897 (1996) 16. J. Habasaki, K.L. Ngai, Y. Hiwatari, J. Chem. Phys. 121, 925 (2004) 17. http://ambermd.org/. The address was confirmed to be valid on 26th Apr. 2016. The latest version is Amber 14. Amber 14 was developed by D.A. Case, V. Babin, J.T. Berryman, R.M. Betz, Q. Cai, D.S. Cerutti, T.E. Cheatham, III, T.A. Darden, R.E. Duke, H. Gohlke, A.W. Goetz, S. Gusarov, N. Homeyer, P. Janowski, J. Kaus, I. Kolossváry, A. Kovalenko, T.S. Lee, S. LeGrand, T. Luchko, R. Luo, B. Madej, K.M. Merz, F. Paesani, D.R. Roe, A. Roitberg, C. Sagui, R. Salomon-Ferrer, G. Seabra, C.L. Simmerling, W. Smith, J. Swails, R.C. Walker, J. Wang, R.M. Wolf, X. Wu and P.A. Kollman (2014), University of California, San Francisco 18. J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kallamn, D.A. Case, J. Comput. Chem. 25, 1157 (2004) 19. I. Todorov, W. Smith, Phil. Trans. R. Soc. Lond. A 362, 1835 (2004) 20. I. Todorov, W. Smith, Phil. Trans. R. Soc. Lond. A 2, 161 (2004) 21. CPMD is available at http://www.cpmd.org/. The address was confirmed to be valid in Oct. 2015 22. S. Plimpton, J. Comp. Phys. 117, 1 (1995). See http://lammps.sandia.gov/, the address was confirmed to be valid in Oct. 2015 23. Thesis of Doctor Courses, C. R. Trott, Technische Universität Ilmenau. GPU code of MD using C was introduced there. 24. http://octa.jp/OCTA/whatsOCTA.html. The address was confirmed to be valid on 19th Feb. 2016 25. http://www.modylas.org/. The address was confirmed to be valid on 19th Feb. 2016 26. http://www.ibiblio.org/pub/languages/fortran/ch1-2.html. The address was confirmed to be valid on 19th Feb. 2016 27. http://www.ks.uiuc.edu/Research/vmd/. The address was confirmed to be valid on 19th Feb. 2016