Polar solvents

As an example, consider the development of molecular mechanics parameters for the nitrile group. A natural choice of appropriate fluid is acetonitrile. First of all we repeat what has been done in article¹⁾.

First of all, we try to reproduce the ΔH_vap = 8.009 kcal/mol and d = 0.773 g/cm³. In the paper calculations were made on MDynaMix. We use the Abalone. It uses slightly different algorithms, so we can expect a small deviation in the results.

Put in a periodic box of 216 molecules of acetonitrile. We take into account the corrections to the pressure (197 atm) and energy (-0.52 kJ / mol = -0.12 kcal / mol), for correction of truncated dispersed interaction. Accordingly, we put the pressure in the barostat at 198 atm, and the energy correction will be added to the heat of vaporization. We choose the same integrator, which was used in the article - MTS VV with 10 short steps in a long one, and a range of “short forces” at 5 Å. Andersen thermostat; a long step 2 fs; simulation time is 200 ps for equilibration and 100 ps for average collection. The average density of the model is equal to 0.779 g/cm³, enthalpy of vaporization is:

ΔH_vap = E_gas - E_liquid + RT

where

R = 8.314472 / 1000 / 4.184 = 0.0019872 kcal/mol K

T = 25ºC + 273.15 = 298.15 K

RT = 0.0019872 * 298.15 = 0.5924 kcal/mol

The average potential in the gas phase was evaluated by simulation (2.5 ns) of a single molecule. PEFRL integrator with a step of 0.5 fs. Andersen thermostat.

E_gas = -22.35 kcal/mol

E_liquid = -6446.8 / 216 -0.12 = -29.97 kcal/mol

ΔH_vap = 8.21 kcal/mol

If intramolecular energy does not vary markedly during the transition from liquid to gas, which is true for most liquids (except water), the heat of vaporization can be estimated by adding RT to the intermolecular energy.

ΔH_vap = -E_{Intermolecular} + RT = -(-1631.7/216 - 0.12) + 0.59 = 8.27 kcal/mol

As you can see, the simplified formula gives almost the same result.

In our calculation both density and heat of vaporization were similar to mentioned in the article. If we have a larger number of molecules, then the correspondence would have been even better.

In order to better feel the situation look how the results depend on the method of simulation. First of all, carry out accurate simulations in order to be able to compare results. We can use any molecular mechanical integrator with a small step size. But it is better to use higher-order integrators, with a step with the middle step size. Fourth-order integrators and PEFRL VEFRL give very good trajectories for the 0.5 fs steps.

The first point shows the density of our model of acetonitrile in the calculation method PEFRL at a step of 0.5 fs, Andersen thermostat. The second point corresponds to the same conditions, but the hydrogen atoms have been given a mass of 10 а.е. In full accordance with the theory (thermodynamic properties are independent of the mass distribution in the system) we have obtained the same result - 0.778 g/cm³. The yellow graph shows the integration step density dependence using the four-order rRESPA integrator. Up until the 10 fs step size density is a constant, but has a small bias to overestimate the density. At the large steps step density is wrong. And with more than 14 steps fs dynamics becomes unstable and the model falls apart. After setting the increased mass of the hydrogen a large step simulation can be carried out (blue graph). Theory suggests that the Monte Carlo simulation should not be depending on the step size chosen. On the brown graph shows the results of simulations using Hybrid Monte Carlo. That is, short dynamic simulation used as a Monte Carlo step (in this case, 100 steps of rRESPA). Indeed, the results do not depend on the magnitude of the dynamic step and correspond to the best estimates by PEFRL. It may seem that these properties make the Monte Carlo method the most effective of the above. However, in the Monte Carlo simulation part of the trajectories are discarded.

Another reason which may decrease the effectiveness of the method lies in the fact that at each Monte Carlo step the velocity of the atoms are updated. This may reduce the diffusion rate and hence the effectiveness of the simulation. Due to these circumstances, molecular dynamics and Monte Carlo methods have a similar efficiency.

We will carry out the optimization of parameters as described in the article, but the charges are developed in the basis aug-cc-pVDZ with the MP2, VDWSCL = 1.1 and with a fixed dipole moment. The results are charges on the nitrogen -0.5236e, on the nitrile carbon 0.4454e, on the methyl carbon -0.2098e, on the hydrogen 0.096e. The density and heat of vaporization have not changed:

d = 0.781 g/cm₃

ΔH_vap = 8.32 kcal/mol

This is an interesting fact. Although the dipole moment of the molecule is not changed, the new charge distribution differs greatly.

We calculated the charge distribution in such a way that the dipole moment of the molecule corresponding to its dipole moment in the gas phase. But we are mainly interested in the liquid state. And in a liquid phase dipole moment is higher than in the gas. According to the literature estimates, it increased by about 15%. In our previous model this effect is not taken into account. The usual practice is to increase the partial charges to simulate the liquid conditions. If all the partial charges increased only by 15% the heat of vaporization and the density sharply increases. But we can compensate this by changing the parameters of the Lennard-Jones potential. If we assume that the methyl group is the same, then we have 4 parameters of nitrile group, which we can change. The number of properties we try to keep two only. It may be that many combinations of these four parameters are suitable. For example, take Lennard-Jones parameters for carbon and nitrogen are the same and will be varying them seeking to a good matching of acetonitrile density and heat of evaporation. Indeed, we can easily find the suitable pair of parameters:

d = 0.768 g/cm³

ΔH_vap = 8.27 kcal/mol

Thus, we have three models of acetonitrile with different parameterizations, but reproducing the experimental values of density and heat of vaporization. Does it mean that such a parameterization completely vague and meaningless? No. All three models can be used in real-world simulations, and they will give better results than the models that were parameterized without reference to the experimental properties of liquids. The fact that this binding is not unique of course somewhat disappointing. But on the other hand, it leaves us room to match our scheme with other molecules, force fields or specific tasks.