Crystalline Silicon (Parrinello-Rahman NPT)

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.

The DFT calculation

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

Molecular Dynamics Parameters

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.[12], 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

Results After 2000 Time-Steps

[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.

Figure 4.4: Fluctuations in the cell angles and lengths during the first 2000 steps of the MD run described in the text.

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.[12] for an example.