Skip to content

The canonical ensemble at constant pressure

Exercise 7: Running NPT Simulations with FHI-aims and I-Pi

Time estimation: \(\sim\)6h

In Part 4 of this tutorial you learned how to perform ab initio MD simulations in the canonical ensemble at constant temperature and volume. However, often experiments are conducted at constant temperature and pressure. To simulate these conditions in theory, the NPT ensemble can be used. The quantities particle number N, pressure P, and temperature T are kept constant during the simulation.

The relevant probability distribution for this ensemble in the isotropic case is usually given by2 $$ \rho_\text{NPT} \propto \text{exp}\left(-\beta\left[K(\vec{p}) + U(\vec{R},V)+P^\text{ext}V\right]\right)$$ where \(K\) is the kinetic energy, \(U\) the potential energy, P the pressure, and V the volume. It should be noted, however, that this definition is not unique1.

As mentioned in Part 4 of this tutorial the pressure can be controlled similarly to the temperature by the introduction of a barostat. The idea behind the barostat is that the physical system can be embedded in an extended system, which exerts the pressure \(P^{\text{ext}}\) such that the internal pressure \(P^{\text{int}}\) is balanced by that. The instantaneous internal pressure \(P^{\text{int}}\) can be calculated from the observables of an ab initio molecular dynamics (AIMD) simulation using the following formula:

\[ P^{\text{int}}_{\alpha\beta} = -\frac{1}{\text{det}(\mathbf{h})} \sum_{\gamma=1}^{3} \sigma_{\alpha\gamma} h_{\beta\gamma} + P^{\text{K}}_{\alpha\beta}, \]

where \(\sigma_{\alpha\gamma}=\frac{\partial E\left[\{\Psi\},\mathbf{R}\right]}{\partial h_{\alpha\gamma}}\) is the stress tensor, \(\mathbf{h}\) is the matrix of lattice vectors, and \(P^{\text{K}}\) is the kinetic contribution to the pressure tensor. However it should be noted that this is not the only

In the following you will learn how to perform an AIMD simulation with constant pressure and temperature using FHI-aims and i-PI.

NPT simulation for copper chloride

Warning

The simulation of CuCl requires computational resources beyond a laptop to obtain statistically converged results in an acceptable time. For your analysis, you can find the trajectories in the solution folder of exercise 7. Please also be aware that we started the simulation from a configuration with very little pressure and that higher pressure simulations can take considerably longer to thermalise.

Copper chloride CuCl has a zincblende crystal structure, thus, the following part will deal with periodic boundary conditions. For our simulation, we will use a 2x2x2 supercell (64 atoms) of the conventional unit cell:

Two notes here:

  1. The above supercell geometry is build from a unit cell that was optimized (lattice parameters and atom positions) with FHI-aims.
  2. The supercell size is another critical numerical parameter, where convergence needs to be tested. For our purposes the 2x2x2 supercell is converged.

For the numerical runtime choices we use the following files for FHI-aims (control.in) and i-PI (init.xml):

Let us have a closer look at the specialities in those input files. As noted in the introduction we need to compute the stress tensor, so i-PI can compute the pressure. In the control.in you additionally need to set:

compute_forces                     .true.
compute_analytical_stress          .true.

For the input.xml file you will need three additional things. First, you have to specify the lattice vectors by setting:

<cell mode="manual" units="angstrom">[10.4366361069840003, 0.0, 0.0, 0.0, 10.4366361069840003, 0.0, 0.0, 0.0, 10.4366361069840003]</cell>

Second, when simulating an NPT ensemble we additionally also have to provide the external pressure:

<ensemble>
    <pressure units='bar'> 0 </pressure>
    <temperature units='kelvin'> 300 </temperature>
</ensemble>

For our simulation, we will use no external pressure, thus, it is set to \(0\).

Third, we need to tell i-PI that we like to run a NPT ensemble. This can be specified in the dynamics section of the input.xml file by setting:

<dynamics mode="npt">
    <barostat mode='isotropic'>
        <thermostat mode='svr'>
            <tau units='femtosecond'> 50 </tau>
        </thermostat>
        <tau units='femtosecond'> 50 </tau>
    </barostat>
    <thermostat mode='svr'>
        <tau units='femtosecond'> 10 </tau>
    </thermostat>
</dynamics>

There are different ways of defining a barostat. A barostat is the equivalent of a thermostat for the pressure: It tries to keep the pressure constant on average. The theory for the barostat defined in the above <barostat> section is given in the paper by G. Bussi, T. Zykova-Timan, and M. Parrinello3.

Now, we are ready to start our simulation. If you do not know how to start the AIMD simulation with i-PI, please read Part 1 of this tutorial. We create four separate folders for each trajectories and choose a maximum of 5000 time steps.

We have prepared a small python script named alltraj.py for you to analyse your trajectories. You find the script in the folder molecular-dynamics-with-i-pi/Tutorial/Exercise-7/solutions/CuCl/, simply copy it to the parent folder of your trajectories and execute it by typing:

python alltraj.py
The script will ask you for the folders in which you have stored your trajectories and up-to which time step you want to plot the results, i.e pressure, temperature, and volume for each trajectory. Additionally the script will also provide the results averaged over
all trajectories.

An example submission script sub.sge for a cluster (using slurm batch system) together with four trajectories in larger unit cells with 5000 time steps, is added to the solution folder, too. The results can be found in the solution folder of exercise 7:

molecular-dynamics-with-i-pi/Tutorial/Exercise-7/solutions/CuCl

This folder will have a subfolder for each trajectory, T1, T2, T3, and T4. In each subfolder, the simulation.out file has the relevant information about the pressure per time step (in column 10). While you are waiting for your trajectories to finish you can use the alltraj.py to analyse the provided simulations. At which point would you consider the system thermalised, can you explain why it thermalised so fast/slow ?

Answer

Since we started from a geometry with almost no initial pressure which is close to our target of zero bar we find that the system thermalises relative quickly. This can best be seen by plotting the pressure and temperature for the first 1000 femto seconds which is shown in the figure below for trajectory T1:

CuCl Pressure Trajectory one CuCl Temperature Trajectory one

The following figure shows the histogram of the computed pressures summed up over all trajectories. However, for each trajectory, we removed the first 2000 steps due to the need of thermalization of the system. To reproduce the figure you can simply execute the script alltraj.py in the folder molecular-dynamics-with-i-pi/Tutorial/Exercise-7/solutions/CuCl/ by typing:

python alltraj.py

Indeed, we find that on average the pressure is zero!

CuCl Pressure


  1. P. Attard, J. Chem. Phys. 103, 9884 (1995). 

  2. [ M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids Oxford University Press, Oxford, 1987.] () 

  3. G. Bussi, T. Zykova-Timan, M. Parrinello, J. Chem. Phys.130, 074101 (2009).