The behavior of hydrogen at extreme pressures and moderate temperatures has relevance to the phase partitioning of giant planetary interiors, as well as the design of systems for inertial confinement fusion. It is challenging to compute hydrogen's equation of state under these conditions, as it is sensitive to the electron dynamics of the system. For example, at higher temperatures, electrons become more diffuse, which cause greater Pauli repulsion and higher pressures. Also, as the nuclei are pressed closer together, hydrogen may attain a metallic character.
We have used eFF to compute the equation of state of bulk hydrogen, and have found it to be consistent with the best theory available over a wide range of temperatures and pressures (J. T. Su and W. A. Goddard III, Phys. Rev. Lett. 99, 185003 (2007)).
In these calculations, we applied periodic boundary conditions to a rectangular unit cell containing 64 hydrogen atoms. For small unit cells, it is tricky to compute the total contribution from long range interactions such as electrostatics, since each charge may see multiple replicas of itself and other charges. The series is convergent, but it is time-consuming to sum these replicas in real space. Instead, eFF uses a modified Ewald algorithm to compute electrostatic interactions, where a portion of the electrostatic energy sum is computed in reciprocal space. The Pauli potential is much shorter ranged, and its periodic sum is computed using a minimum image convention.
To apply periodic boundries, set the option periodic = true, and specify the rectangular boundaries with the parameters x_bound, y_bound, and z_bound. The current Ewald implementation is rather slow, so bear in mind that computations with this option can be as much as ten times slower than the equivalent nonperiodic calculation.
As an example, we compute the pressure of hydrogen at 15,300 K and = 2 bohr, where the density is 1 . The operation proceeds in two stages. First, we optimize hydrogen at = 2 bohr and 0 K by building a lattice of hydrogen molecules, minimizing them, annealing them at a moderate temperature, and minimizing them again. Next, we heat the hydrogen up to 15,300 K and compute the average pressure and temperature over the course of the dynamics run.
To fill the unit cell with hydrogen molecules, we use the script hydrogen_cfg [nx] [ny] [nz] [rs], which places hydrogen atoms in a unit cell, arraying them on a regular rectangular lattice, with alternating layers displaced in a way that causes the atoms to pair into molecules upon minimization. We then optimize the geometry:
% h2bulk_cfg 4 4 4 2.0 > h2bulk_min.cfg % eff h2bulk_min.cfg h2bulk_min.cfg: @params x_bound = 0.000000 12.895936 y_bound = 0.000000 12.895936 z_bound = 0.000000 12.895936 periodic = true calc = minimize ...
We anneal the hydrogen box to find a lower-energy local minimum, since the hydrogens are in an unstable configuration. We do this using NVT dynamics with a Nose thermostat (thermostat = nose_hoover), ramping the temperature linearly down from start_temperature = 3,000 K to end_temperature = 0 K.
% cp h2bulk_min.cfg.restart h2bulk_anneal.cfg % [add NVT parameters to h2bulk_anneal.cfg] % eff h2bulk_anneal.cfg h2_bulk_anneal.cfg: @params x_bound = 0.000000 12.895936 y_bound = 0.000000 12.895936 z_bound = 0.000000 12.895936 periodic = true calc = dynamics thermostat = nose_hoover start_temperature = 3000.0 end_temperature = 0.0 dt = 0.02 num_steps = 10000 print_every = 100 ...
We then reminimize the annealed structure:
% cp h2bulk_anneal.cfg.restart h2bulk_min2.cfg % [add minimization parameters to h2bulk_min2.cfg] % eff h2bulk_min2.cfg h2bulk_min2.cfg: @params x_bound = 0.000000 12.895936 y_bound = 0.000000 12.895936 z_bound = 0.000000 12.895936 periodic = true calc = minimize ...
With a reasonable minimum in hand, we perform NVE dynamics on the system to obtain the equation of state data point. We initialize the temperature to 30,600 K - since the energy is at a minimum, the virial theorem suggests that after equilibration, the average temperature will be close to 15,300 K.
% cp h2bulk_min2.cfg h2bulk_EOS.cfg % [add NVE parameters to h2bulk_EOS.cfg] % eff h2bulk_EOS.cfg h2_bulk_EOS.cfg: @params x_bound = 0.000000 12.895936 y_bound = 0.000000 12.895936 z_bound = 0.000000 12.895936 periodic = true calc = dynamics start_temperature = 30600.0 dt = 0.02 print_every = 50 num_steps = 100000 ...
The program automatically computes instantaneous temperature and pressure over the dynamics run using virial expressions, e.g. and , where we sum over all nuclear and electronic degrees of freedom and take .
The temperature and pressure are reported in lines containing the tag [dynamics_iter]:
% cat h2_bulk_EOS.eff | grep dynamics_ ... [dynamics_header] ... measured_temp p_ke_rigid p_ke_flexible ... [dynamics_iter] ... 25701.454374 0.000000e+00 0.000000e+00 ... [dynamics_iter] ... 17807.946509 4.862479e+10 4.995476e+10 ... [dynamics_iter] ... 16673.384279 3.809880e+10 4.622077e+10 ... [dynamics_iter] ... 16017.964512 3.711733e+10 4.478320e+10 ... [dynamics_iter] ... 10704.564113 2.670103e+10 2.998048e+10 ... ...
The pressure is reported as four components, p_ke_rigid, p_ke_flexible, p_pe_rigid, and p_pe_flexible. The kinetic energy (KE) portion is the term from the ideal gas law, and is proportional to temperature. To understand its origin, consider what happens when a system is squeezed. The particle velocities appear to increase relative to the smaller size of the system, which in turn increases the pressure. The potential energy (PE) portion is the term that accounts for intermolecular interactions, and is dominated by the Pauli repulsion between molecules.
Rigid pressures assume that the electron size stays fixed when the system is squeezed, while flexible pressures assume that the electron size shrinks uniformly along with the system. Which pressure is correct depends on the extent to which electrons are associated with nuclei. In a gas of low density helium atoms, the rigid pressure is correct, while in a uniform electron gas or metal, the flexible pressure is correct. In general, the true pressure will be in between these two extremes. Usually the flexible and rigid pressures are similar, and in this case we select the rigid pressure to compute our EOS data point.
We throw out the first 0.5 ps of dynamics, since the system is equilibrating to the final temperature during that time. Using a spreadsheet to average together temperatures and pressures from the last 1.5 ps, and summing together kinetic energy and potential energy contributions to the pressure, we obtain = 12,986 625 K, = 55.6 6.6 GPa, and = 55.6 3.7 GPa. The error bars may be reduced further by running a longer simulation.
Figure 2.12 shows an EOS curve obtained by repeating the above calculation for a variety of temperatures. There is excellent agreement between the eFF data and those obtained from other higher level theory.