STRUCTURE OF MATTER AND QUANTUM CHEMISTRY Molecular Dynamics Simulation of Diisopropyl Ether Using Various Interatomic Potentials

A comparative assessment of the accuracy of determining the density and viscosity has been carried out for diisopropyl ether using the method of classical molecular dynamics using three potentials. The accuracy of determining the viscosity coefficients when using equilibrium and nonequilibrium calculation methods was also investigated.

Diisopropyl ether is a simple aliphatic ether with the formula C 6 H 14 O, a colorless liquid that is soluble in organic solvents. It is used as an organic solvent in chemical technology in the production of dyes and resins, as well as for the isolation and purification of various substances by extraction [1][2][3]. The addition of diisopropyl ether to gasoline, similarly to the addition of methyl tert-butyl ether, makes it possible to reduce the emissions of soot particles and carbon monoxide during the operation of internal combustion engines [4][5][6][7].
Diisopropyl ether is of practical interest not only as an industrial solvent and fuel additive, but also as a component of electrochemical current sources and systems for extracting rare metals from solutions. Ethers can be a key component of ion-selective barriers, which are a "liquid membrane" [8], where a layer of an organic liquid immiscible with water separates two layers of an aqueous solution. The selective permeability of such a layer with respect to different ions can be used both to isolate rare elements such as lithium or rubidium from aqueous solutions, and to create new promising electrochemical current sources such as Red/Ox flow elements [9,10].
Ethers: crown ethers and cryptands are among the most interesting components that can provide the selectivity of a "liquid membrane." Moreover, while there is a fairly wide set of experimental data for diisopropyl ether in the literature, for crown ethers and cryptands they are much more scarce. This makes diisopropyl ether a convenient object for developing molecular dynamics modeling methods, which in the future should be suitable for a wider class of ethers, including crown ethers and cryptands.
When describing the operation of an ion-selective barrier such as a liquid membrane, the description of the mobility of processes in an ether medium and solutions containing it is of the greatest importance. Moreover, the quality of the description of the transfer processes correlates well with the accuracy of the description of the transport properties of the system: viscosity and diffusion.
Although diisopropyl ether is a widely used and readily available substance, the study of the properties of its mixtures with other organic substances is still an urgent task [6,7,11]. The variety of mixtures of organic substances, including diisopropyl ether, makes it expedient to study them using computer chemistry methods, such as the molecular dynamics (MD) method, which are widely used to study the structure and properties of liquids [12][13][14][15]. Examples of calculating the dynamic properties of aqueous solutions are given in [12,[16][17][18].
The physical basis of the methods of molecular dynamics are interatomic potentials, that is mathematical functions for calculating the potential energy of a system of atoms with given positions in space. Therefore, the accuracy of modeling when using the MD method is determined by how correctly the used interatomic potentials reproduce the force interaction between atoms and molecules of the substance under study. Verification of a wide variety of interatomic potentials based on experimental or quantum chemical data is an important task of modern MD [19][20][21][22].

STRUCTURE OF MATTER AND QUANTUM CHEMISTRY
It should be noted that the works devoted to the molecular dynamics modeling of diisopropyl ether and its mixtures with other substances are not numerous [23,24]. In most of them, the authors conducted research using only one force field. Thus, the relevance of this work lies in assessing the applicability of the three classical potentials GAFF (AMBER), OPLS-AA, and CHARMM in terms of the accuracy of reproduction by molecular dynamics methods of the transport properties of diisopropyl ether in a wide range of temperatures and pressures.
In Section 2 of this work, the features of the used interatomic interaction potentials are considered (Section 2.1), the main parameters and characteristics of molecular dynamics modeling are given (Section 2.2), and the methods used to calculate the viscosity are described (Section 2.3). Section 3 presents the results of estimating the density (Section 3.1) and compressibility (Section 3.2) of diisopropyl ether depending on pressure and temperature (Section 3.3), the accuracy of determining the viscosity coefficients when using various methods for their evaluation is considered (Section 3.4), a comparative assessment of the accuracy of determining the transport properties of diisopropyl ether using three classical potentials (Section 3.5) is performed.

Potentials of interatomic interaction
For modeling diisopropyl ether, three interatomic potentials were considered for applicability: General Assisted Model Building with Energy Refinement Force Field (GAFF) [25], Optimized Potentials for Liquid Simulations All-Atom (OPLS-AA) [26] with charge correction, and Chemistry at HARvard Macromolecular Mechanics (CHARMM) [27], version 36.
The parameterization of interactions in the GAFF potential was formed using the Antechamber program. The calculation of partial charges on atoms was carried out according to [28], the parametrization of interactions was carried out according to [29,30]. The parameters for OPLS-AA/CM1A were generated on the LigParGen server [31]. For partial charges, the correction 1.14*CM1A was used [32,33]. The interaction constants in CHARMM36 were generated in the CHARM-GUI [34][35][36].
All three potentials include covalent and nonvalent interactions between atoms: For GAFF and OPLS-AA/CM1A, covalent interactions are described by harmonic vibrations of covalent bonds, angles between three atoms, and torsion interactions: These terms are defined as follows: are the energy constants, -l 0 , θ 0 are the equilibrium values of the bond length and the angle between atoms, respectively, and ϕ is the torsion angle between the corresponding planes drawn through the triples of atoms.
For the CHARMM36 potential, the formula for the energy of covalent interactions E valent also contains the terms E bond , E angle , and E dihedral , which are described by the same functional dependencies. In addition, it contains two more terms, that is, E UB and E improper . The term E UB is called the Urey-Bradley sum. In it, summation occurs over all chains of three bonded atoms A-B-C. The E improper parameter describes the out-of-plane bending energy of atoms, and the summation occurs over a set of four atoms that are not connected in series. These terms have the following form: where k u , k imp are the energy constants, ψψ 0 is the dihedral angle formed by four atoms not connected in series, and s is the distance between atoms A-C, in the chain of atoms A-B-C.
For all three potentials, nonvalent interactions of atoms are described by the Lennard-Jones and electrostatic potentials: where r ij is the distance between interacting atoms, ε and σ are the Lennard-Jones parameters, and q i and q j are the atomic charges.
For the GAFF and CHARMM36 potentials, the parameters ε and σ are determined using the Lorentz-Berthelot combining rule [37,38], while the geometric mean combining rule is used for the OPLS-AA/CM1A potential, [39,40].
For atoms within the same molecule separated by one or two bonds, nonvalence interactions are not taken into account. Nonvalence interactions are taken into account in a special way for atoms separated by three bonds: in OPLS-AA/CM1A for these atoms, the Lennard-Jones interaction and the electrostatic interaction are scaled with a factor of 0.5; in GAFF, elec- trostatic interactions and Lennard-Jones interactions are scaled with coefficients of 0.833 and 0.5, respectively; in CHARMM36, the electrostatic interaction is not scaled, and the Lennard-Jones interaction is determined based on a special set of parameters.

Methods of Molecular Dynamics
The simulation was carried out using the GROMACS software package [41]. The simulation was carried out for a cell with 3375 molecules. The initial configuration was set in the form of a cubic lattice 15 × 15 × 15 molecules in size. Then, the cell was compressed to a density corresponding to the experimental one for diisopropyl ether at the given temperature and pressure. After compression, the system was brought to the desired temperature in the NVT ensemble and then brought to the desired pressure in the NPT ensemble. Both processes were simulated over an interval of 200 ps.
To maintain the NVT ensemble, a modified Berendsen thermostat [42] was used, while for the NPT ensemble a Berendsen barostat was used [43]. After equilibrium was established in the NPT ensemble, the equilibrium density was calculated. In all calculations, the time integration step was 1 fs. Periodic boundary conditions were used to exclude edge effects.
In the simulation, the Lennard-Jones and electrostatic interaction potentials were trimmed. In this case, corrections to the energy and pressure were taken into account, which compensated for the cutoff of potentials [44]. To assess the effect of the cutoff radius on the simulation results, calculations were performed for two values of this parameter: 1 and 1.2 nm. It was found that with an increased cutoff radius, the accuracy of modeling the dynamic and other properties of diisopropyl ether does not increase; therefore, the results are presented in the work only for 1 nm.
The long-range part of the Coulomb potential was calculated by the Ewald method (SPME) [45].

Methods for Calculating Viscosity in MD
The Green-Kubo method [46,47] was used as the main calculation method: the viscosity value is determined from the relation where V is the volume of the computational cell, k is the Boltzmann constant, T is the absolute temperature of the system, and C σ is the autocorrelation function of the off-diagonal elements of the stress tensor Angle brackets mean ensemble averaging and three mutually perpendicular planes xOy, xOz, yOz. The stress tensor is calculated by the formula: where -β is the component of the force acting on the ith particle, N is the number of atoms in the system, and N' is the number of atoms in the system and the nearest image in the case of periodic boundary conditions.
To calculate the viscosity coefficients by this method, simulation was carried out in the NVT ensemble for 15 ns. Then, 150 statistically independent autocorrelation functions C σ were calculated from the obtained trajectory, which were then averaged, and the viscosity was calculated from the obtained function.
Two additional viscosity evaluation methods were tested for applicability: the Stokes-Einstein method and the nonequilibrium method.
The Stokes-Einstein method is based on the Einstein relation: where R is the effective radius of the molecule, D is the self-diffusion coefficient, k is the Boltzmann constant, and T is the absolute temperature of the system. The average molecular gyration radius was taken as the effective radius R.
The self-diffusion coefficient D was calculated using the Einstein-Smoluchowski formula, which describes the dependence of the root-mean-square displacement (RMSD) of molecules at long times: where r is the radius vector of the particle coordinate at time t, r 0 is the radius vector of the particle coordinate at the moment corresponding to the origin of time.
The RMSD was calculated in the NVT ensemble within 1 ns. In this case, several RMSD were calculated from one trajectory: the initial reference points were taken along the entire trajectory with a period of 10 ps, the RMSD was calculated for them, and then averaging was carried out.
The nonequilibrium method for calculating viscosity belongs to the methods of nonequilibrium molecular dynamics. Atoms of substance are placed in an external force field of the form [48].
In such a system, viscosity is estimated as the response of a substance to shear stress. Since the external force does work on the displacement of atoms, the thermodynamic state of the system is disturbed. This can lead to a significant increase in pressure [48]. Minimization of this phenomenon is achieved by selecting the optimal values of the coefficient A and the size of the computational cell L x , L y , L z . By enu- No. 6 2023 meration of various parameters, it was found that for this system with the GAFF potential, the optimal values are as follows: A = 0.005 nm/ps 2 , L x = 5.73 nm, L y = 6.02 nm, L z = 23.4 nm at 3500 molecules. With these parameters, the viscosity coefficients were calculated. The NVT ensemble was used and the simulation was carried out for 5 ns.

Pressure Dependence of Diisopropyl Ether Density
At the first stage of the work, the density at a temperature of 303.15 K was calculated for various pressures using the considered potentials. Based on the data we obtained, a graph of the dependence of density on pressure was plotted (Fig. 1). The experimental density values are taken from [49].

Compressibility of Diisopropyl Ether
Calculation of the compressibility factor β of diisopropyl ether in the studied range of pressures and temperatures was carried out according to the formula: where ρ is the density of the system and P is the pressure in the system.
Since the dependence of pressure on density is close to linear in the pressure range under study, the derivative was calculated as the slope of the approximating straight line.
To estimate the compressibility factor, we used the density values obtained by the MD method using the CHARMM36 potential, since this potential provided the highest accuracy. The calculation results are shown in Table 1.

The Temperature Dependence of Diisopropyl Ether Density
The equilibrium density was calculated at a pressure of 0.1 MPa for temperatures in the range from 243.15 to 333.15 K using the considered potentials. Figure 2 shows the obtained values together with the experimental data [49].

Viscosity Estimates Using Different Methods
The viscosity coefficients were calculated using three different methods in the temperature range from 273.15 to 333.15 K and pressure 0.1 MPa using the GAFF potential.   The dependence of the viscosity coefficient η on temperature is well described by the Arrhenius equation: where E a is the activation energy, R is the universal gas constant, T is the absolute temperature of the system, and η 0 is the preexponential factor. Accordingly, the logarithmic dependence of viscosity on the reciprocal temperature is a linear function. Figure 3 shows plots of viscosity versus reciprocal temperature, on a logarithmic scale.

Calculation of Viscosity Coefficients Using Various Potentials
For the same temperatures and pressures, as in the previous Section, the viscosity coefficients were calculated using the GAFF, OPLS-AA/CM1A, and CHARMM36 potentials. The Green-Kubo method was used as the most common and convenient method for calculating viscosity. The calculation results are shown in Fig. 4.

RESULTS AND DISCUSSION
Using the molecular dynamics method using the GAFF, OPLS-AA/CM1A, and CHARMM36 potentials, some properties of diisopropyl ether were calculated for various temperatures and pressures.
1. When assessing the dependence of the density of diisopropyl ether on pressure (Fig. 1), the obtained values turned out to be somewhat higher than the experimental ones. The GAFF and CHARMM36 potentials give similar values: for these, the deviation does not exceed 20 kg/m 3 (3%). For OPLS-AA the deviation does not exceed 30 kg/m 3 (5%).
2. Evaluation of the compressibility of diisopropyl ether at a temperature of 303.15 K (Table 1) showed good agreement with the experimental data. When using the CHARMM36 potential, the deviation of the results of molecular modeling from experiment does not exceed 12%.
3. Calculation of the density of diisopropyl ether was carried out at ten points in the temperature range from 243.15 to 333.15 K using all three potentials. The calculation results and experimental data are shown in Fig. 2. The best calculation accuracy is provided by the GAFF and CHARMM36 potentials, which gave values close to each other. For them, the deviation from experiment does not exceed 25 kg/m 3 (3.3%). The OPLS-AA/CM1A potential provides the least accuracy: the error is no more than 40 kg/m 3 (5.1%).
4. The calculation of the viscosity coefficient was carried out using three different methods. The need to use three methods was due to the fact that the initial calculations of the viscosity coefficient of diisopropyl ether using the Green-Kubo method showed a rather large discrepancy with the experimental data. Therefore, it became necessary to evaluate this parameter using other methods to make sure that the error is caused by insufficient accuracy of the mathematical model and not by the viscosity estimation method. The results of calculating the viscosity using 3 methods and the experimental data are shown in Fig. 3.
One can see that all three methods gave similar values. This confirms that the viscosity estimation errors are not associated with an error in the methods of its calculation, but with the inaccuracy of the mathematical model caused by the inaccurate correspondence of the  potentials used to the real force fields that act inside the substance. 5. Comparison of the results of viscosity calculation using three potentials (Fig. 4) shows that the CHARMM36 potential provides the highest viscosity modeling accuracy. The viscosity values calculated using this potential differ from the experiment by 20-70%, and the higher the temperature is, the higher the accuracy is. The viscosity estimation errors using GAFF and OPLS-AA/CM1A potentials are in the range of 50-140%. Thus, we can recommend the use of the CHARMM36 potential for studying the transport properties of diisopropyl ether and its mixtures with other ethers if high accuracy of calculations is not required.
Thus, the CHARMM36, GAFF, and OPLS-AA/CM1A potentials reproduce the dependence of the density and compressibility of diisopropyl ether on temperature and pressure well.
Among the investigated potentials, the CHARMM36 force field gives the best description of the viscosity of diisopropyl ether and its dependence on temperature. This suggests that this potential should also give a satisfactory description of other transfer processes in diisopropyl ether, as well as in mixtures of diisopropyl ether with other ethers, among which solutions of crown ethers in diisopropyl ether are of the greatest practical interest in connection with the task of creating ion-selective barriers of the liquid membrane type.
However, it should be noted that all three potentials give a noticeable discrepancy between the experimental values of the viscosity of diisopropyl ether and the results of molecular dynamics simulations. One of the reasons for this discrepancy may be that the parameters of nonvalent interactions on atoms in an organic molecule are universal for an atom with each type of bond and do not fully take into account the inductive effect of various substituent groups. It is possible that, while maintaining the structure of the considered potentials, redefining the nonvalent interactions of atoms from those given in the standard form to those calculated individually for the molecule under study can help improve the agreement between the results of molecular dynamics simulation and the actual behavior of the substance in the experiment. However, this issue is a task for further research.