PES for air molecules

Potential energy surface of interaction of two diatomic molecules for air flows simulation

1. DFT data

For each of the three types of interaction, DFT calculations were performed for the following mutual orientations of the molecules (see figure):

with equilibrium interatomic distances b. Thus, 18810 points were obtained for each pair of molecules (N2‒N2, N2‒O2, O2‒O2), a total of 56430 DFT calculations.

Download DFT points: N2-N2.txt, N2-O2.txt, O2-O2.txt



2. Two-particle interaction model

An attempt was made to describe the potential energy surface (PES) obtained in the DFT calculations using pair interactions between the atoms of different molecules. This was done because this case of the interaction model is very simple and very convenient for trajectory MD calculations.

It was obtained that it is impossible to describe the interaction energy of the studied diatomic molecules using any two-particle interaction model with a mean-square relative discrepancy e_rel≲ 6% (even in the considered simplified case of rigid molecules). In this regard, below we tried to describe the obtained DFT points using the multi-particle interaction model.


Obtained mean-square relative discrepancy e_rel for two-particle fitting models

∆r, A N2‒N2 N2‒O2 O2‒O2
0.1 7.86%
Download
5.8%
Download
5.74%
Download
0.025 7.61%
Download
5.5%
Download
5.58%
Download


3. Multi-particle interaction model

3.1. Model description

The potential was divided into the short-range (SR) and long-range (LR) parts: U = U_SR + U_LR. The SR part was U_SR = ∑ U_M(r_i), where U_M(r) = D*(e^(-2*α*(r-a)) - 2*e^(-α*(r-a))) ‒ Morse potential, α, a, D ‒ constants, r_i ‒ distances between atoms of different molecules. The LR part was written as a combination of trigonometric functions of the orientation angles: U_LR = A(R)*(∑ L_ijkl * F_i(γ_0) * F_j(γ_1) * F_k(γ_2) * G_l(R)), where L_ijkl ‒ constants, F_i(γ_j) ‒ i-th element of the set B_j = {1, cos(γ_j), sin(γ_j), cos(2*γ_j), sin(2*γ_j), } for j = 0,1,2, G_i(R) ‒ i-th element of the set B_3 = {1, R^(-1), R^(-2), }, N_j ‒ number of elements in the corresponding set B_j, A(R) = (R/R_0)^β exp(-β*R/R_0), where R_0, β ‒ constants.


Main characteristics of the proposed fitting models

Model 1 Model 2
Feature N2‒N2 N2‒O2 O2‒O2 N2‒N2 N2‒O2 O2‒O2
α 2.5 2 2.5 2 3.5 2.5
β 3 4 4 4 5 5
R_0, A 1.3 1.4 1.4 1.3 1.8 1.6
Number of terms 377 74
e_rel 0.73% 0.68% 0.73% 1.88% 1.56% 1.70%
K_j DownloadDownloadDownloadDownloadDownloadDownload
Init file m1-N2N2.txtm1-N2O2.txtm1-O2O2.txtm2-N2N2.txtm2-N2O2.txtm2-O2O2.txt


3.2. Model code (C++)

The described models and force field computation were implemented as a C++ code:

Before use, it is necessary to initialize the model from the file:

        // Model class
	airpotentials model;

	// Init model from file
	model.LoadModelFromFile("m1-N2N2.txt");

To calculate the energy and forces, it is necessary to call function ComputeEnergy(...) with the coordinates of the atoms and a pointer to the array into which the forces on the atoms will be stored:

        // Atoms coordinates and forces
        double x[4][3] = {...};
	double f[4][3];
	
        // Energy and forces computation
        double U = model.ComputeEnergy(x,&f[0][0]);


Acknowledgments

The research is carried out using the equipment of the shared research facilities of HPC computing resources at Lomonosov Moscow State University and Joint Supercomputer Center of the Russian Academy of Sciences. The research was supported by the Russian Foundation for Basic Research (project 18-31-20025).