Applications of High Performance Computing: Born–Oppenheimer Molecular Dynamics of Complex Formation in Aqueous Solutions

The progress of supercomputer technologies initiated the development of methods of computational chemistry and their applications, particularly molecular dynamic simulations with ab initio potentials. These new methods allow to solve important problems of chemistry and technology. Particularly, solvent extraction and separation techniques are widely used to decrease the amount of radioactive wastes, especially radioactive caesium isotopes present in liquid phases. We demonstrated that the calculated binding constants between the alkali cation and calix[4]arene differ 10 times for Cs and Na ions, that is in good agreement with the experimental value. We report the results of benchmark calculations of our model system composed of 929 atoms described in the density functional theory approximation with the GGA-type functional PBE with the empirical dispersion correction D3 and combined basis of Gaussian functions and plane waves DZVP with Goedecker-Teter-Hutter pseudopotentials. We demonstrate that efficiency of calculations decrease to about half if the amount of nodes is 16 on the Lomonosov-2 supercomputer.


Introduction
The progress of supercomputer technologies initiated the development of methods of computational chemistry and their applications. In this paper we focus on the calculations of the free binding energies of the complex formation between alkali metal cations and calix [4]arene ligand. Calix [4]arene and its derivatives are successfully used in solvent extraction and separation techniques to decrease the amount of radioactive wastes containing caesium isotopes. In the end of 20 th century computational works aiming to study these interactions were based on either classical molecular dynamics simulations [6], or ab initio calculations in the gas phase [1,2]. The first approach suffers from the shortcomings of the classical force fields especially when analyzing interactions of the charged particles like metal cations. The second group of methods can better describe the interaction, although they do not account for the entropic contributions important for the binding process and do not explicitly consider solvent molecules.
In the present study we perform molecular dynamics simulations with the DFT (density functional theory) potentials to study the binding processes of sodium and caesium cations with the calix [4]arene ligand considering explicitly the solvation water box. We also perform benchmark tests to evaluate the computational cost of these simulations and suggest the optimal amount of the computational resources for this type of calculations.

Models and Methods
The model system comprised 291 water molecules, calix [4]arene with one of the four alcohol groups being deprotonated, and the alkali metal cation (Na + or Cs + ), 929 atoms in total. It was preliminaryly equilibrated to the constant volume of 20.8 × 19.4 × 22.3Å 3 (Fig. 1) in classical molecular dynamic (MD) simulations. The Born-Oppenheimer MD was performed in the CP2K program package [8].

Me + calix[4]arene
Me + R Figure 1. The model system (left) and the reaction coordinate between the alkali metal cation and oxygen atom of calix [4]arene (right). Carbon atoms are shown in green, oxygen in red, metal cation in magenta and hydrogen in white The entire system was described in the density functional theory approximation with the GGA-type functional PBE [7] with the empirical dispersion correction D3 [5] and combined basis of Gaussian functions and plane waves DZVP with Goedecker-Teter-Hutter pseudopotentials [3]. The calculations were performed with the biasing harmonic potentials centered at the certain values of the reaction coordinate being the distance between the alkali metal cation and one of the oxygen atoms of calix [4]arene. For the system containing caesium cation 16 simulations with the different binding potentials were chosen, and for the system with the sodium cation -8. Each trajectory was 10 ps with the integration time step of 0.5 fs at T=300 K. The trajectory analysis was performed with the WHAM (weighted histogram analysis method) and UI (umbrella integration) approaches.

Results and Discussion
We begin with the benchmark calculations of our model system on the Lomonosov-2 supercomputer (1696 nodes, single-socket Intel Haswell-EP E5-2697v3, 64Gb RAM, NVidia Tesla K40M) [8] (Tab. 1). It is more important to focus on the strong scaling efficiency of the single MD step rather than single SCF step as MD simulations are performed. If the amount of nodes is 16, the efficiency decreases to about half therefore we did not use more than 12 nodes in our calculations (see Tab. 1). Figure 2 demonstrates the standard Gibbs energy profiles calculated for the sodium and caesium cations. The equilibrium distances between the oxygen atom of the alkali cation and the oxygen atom of calix [4]arene are 2.6Å and 3.5Å for sodium and caesium cations, respectively. It is in line with the ionic radius difference being around 0.7Å. The calculated Gibbs free energy of complex formation is approximately 6 kcal/mol larger in case of caesium cation resulting in approximately 5 × 10 4 difference in dissociation constants, that is in a good agreement with the experimental value of 5 × 10 3 [4].

Conclusion
We demonstrate that modern supercomputer facilities allow calculating free energy profiles of metal cation -ligand binding in aqueous solution. MD simulations with the DFT potentials implemented in the CP2K program package demonstrate good scalability and correct description of the model systems.