In our second example, we will use Langevin Parrinello-Rahman dynamics to investigate structural fluctuations in silicon. While each time-step is computationally cheap, this example will require many thousands of MD steps to obtain useful results.
We will use the following cell file for the MD calculation. This represents the geometry optimised structure converged with respect to plane wave basis size at 250 eV, and with respect to the k-point grid at 2x2x2 k-points. An offset in the Monkhorst-Pack grid appropriate for the diamond structure has been used. An external pressure of approximately atmospheric has been applied.
%block LATTICE_ABC Ang 5.37639 5.37639 5.37639 90.0 90.0 90.0 %endblock LATTICE_ABC %block POSITIONS_FRAC Si 0.000000 0.000000 0.000000 Si 0.000000 0.500000 0.500000 Si 0.500000 0.500000 0.000000 Si 0.500000 0.000000 0.500000 Si 0.750000 0.250000 0.750000 Si 0.250000 0.250000 0.250000 Si 0.250000 0.750000 0.750000 Si 0.750000 0.750000 0.250000 %endblock POSITIONS_FRAC %block SPECIES_POT Si Si_00.usp %endblock SPECIES_POT %block external_pressure MPa 0.1 0.0 0.0 0.1 0.0 0.1 %endblock external_pressure kpoints_mp_grid 2 2 2 kpoints_mp_offset 0.125 0.125 0.125 fix_com=true
Again we follow a similar procedure to that given in section 4.1.2 to identify a suitable time-step of 2 fs. We also need to identify suitable relaxation times for the thermostat and barostat. As discussed in ref., we may wish to compute a bulk memory function using an empirical potential for silicon, as well as studying the phonon density of states, to compute parameters which will most accurately reproduce the bulk dynamics. For this example we will simply use ``best-guess'' parameters which will sacrifice some accuracy in short-term dynamics in favor of increased sampling efficiency.
cut_off_energy = 250 eV fix_occupancy = true task = MOLECULARDYNAMICS md_num_iter = 500000 md_delta_t = 2 fs md_ensemble = NPT md_temperature = 293 K md_ion_t = 7.0 ps md_cell_t = 20.0 ps md_thermostat = LANGEVIN md_barostat = PARRINELLO-RAHMAN fixed_npw = false
[N.B. with long calculations such as this, the unix command
grep -c '<-- E' myrun.md
is very useful for counting how many time-steps have been written to the .md file.]
After 2000 steps the system has begun to settle into its equilibrium fluctuations of temperature and density. Fluctuations in the cell parameters are shown below.
Note that this number of time-steps is insufficient to produce a correctly converged distribution of temperature and volume samples. With many thousands more steps, it is possible to calculate heat capacity, bulk modulus, and expansion coefficient from these distributions. See ref. for an example.