Viscosity: Green-Kubo relation

This tutorial will show how to calculate the viscosity of liquid benzene using the Green-Kubo relation:

\[\eta = \frac{V}{k_{\rm{B}}T} \int_0^\infty <P_{\alpha\beta}(0)P_{\alpha\beta}(t)> \rm{d} t\]

where \(\eta\) is the viscosity, \(V\) is the cell volume, \(k_{\rm{B}}\) is the Boltzmann constant, \(T\) is the temperature, \(P_{\alpha\beta}\) is an off-diagonal component of the pressure tensor, and \(t\) is time. The function \(<P_{\alpha\beta}(0)P_{\alpha\beta}(t)>\) is an autocorrelation function with the average taken over all time origins and off-diagonal components (αβ is yz, xz, or xy).

The pressure tensor equals the stress tensor plus a kinetic contribution. This tutorial was also featured in our video tip of the week series

Step 1: Create a box of benzene

Create a box of 16 benzene molecules with density 0.875 g cm^-3 (experimental room-temperature density of benzene) using the Builder in AMSinput, or download and import a premade xyz file.

Start AMSinput

Edit → Builder

Enter 13.333 on the diagonal to change the box to a cube of 2370.2 ų

Next to Fill box with:, enter 16

Next to copies of: start typing benzene and then click the Benzene (ADF) entry in the list that pops up

Click the Generate Molecules button

Click the Close button at the bottom of the Builder window

The molecules are generated at random positions and orientations, with constraint that all atoms (between different molecules) are at least the specified distance (2.5 Å) apart.

Step 2: Set up the benzene MD simulation

The next step is to set up the details of the simulation.

Switch to the ForceField panel: →

Set Periodicity → Bulk

Set Type → GAFF

Check the Automatic atom typing checkbox

(optional) Edit → Crystal → Map atoms to 0..1 will place the atoms inside the unit cell

At the bottom right, click the Pre-optimize button and wait for the pre-optimization to finish

Set Task → Molecular Dynamics

For the MD simulation, we will set up a room-temperature simulation thermostatted with a Nose-Hoover thermostat. It is important to set Sample frequency to a small number. It will save not only the atomic coordinates but also the pressure tensor at that sample frequency, and the pressure tensor needs to be saved frequently to calculate the integral accurately.

Model → MD or next to Task: Molecular Dynamics

Number of steps: 100000

Time step: 1 fs

Sample frequency: 5 (or smaller)

Initial temperature: 298 K

Click on next to Thermostat

Click on the to add a thermostat

Thermostat: NHC

Temperature(s): 298

Damping constant: 100 fs

We must also tell the AMS driver to calculate the pressure (if you use a barostat, that will happen automatically). To save disk space, we can also turn off most of the things usually saved to the trajectory file, e.g. velocities and molecular bonds, as well as the engine checkpoint files.

Model → MD

Set Checkpoint frequency to 100000

Click the next to Trajectory Options

Check Calculate pressure

Uncheck Write velocities

Uncheck Write charges

Uncheck Write bonds

Uncheck Write molecules

Step 3: Run the simulation

Now we will run your set-up:

Use the File → Run command

When asked to save your input, save it with the name Benzene

Important

Double-check that the pressure is calculated. At the start of the simulation, open the logfile (SCM → Logfile), and make sure that numerical values (and not “n/a”) is printed for the pressure. If there is no pressure, then stop the simulation and go back to Step 2.

Note

Running this calculation takes approximately 20 minutes.

Step 4: Analyze the simulation results

Wait until the calculation has finished, and then open it in AMSmovie.

Start AMSmovie: SCM → Movie in the logfile window

Press the Play button (the triangle pointing right at the left bottom of the AMSmovie window)

The pressure tensor autocorrelation function should only be calculated after the simulation has equilibrated.

In the AMSmovie window:

Use the MD Properties → Potential Energy command

Use the MD Properties → Time command

The simulation seems to equilibrated quite quickly. After frame 4000 (time 20 ps) the potential energy seems to fluctuate around a constant value, indicating that the simulation as likely equilibrated at that point.

Frame 4000 corresponds to 20 ps because the time step was 1 fs and the sample frequency 5: 4000*5*1 fs = 20000 fs = 20 ps.

MD Properties → Autocorrelation function

Change the starting step from 1 to 4000 (this will discard the first 20 ps of the simulation)

Property → Pressure Tensor

Vector elements: 4 5 6 (these are the off-diagonal elements \(P_{yz}\), \(P_{xz}\), and \(P_{xy}\))

Max ACF Step: 3000 (corresponding to a maximum correlation time of 15 ps)

Press the button Generate ACF

Close the Autocorrelation Function window

This brings up several graphs in AMSmovie. The bottom graph contains the integral of the pressure tensor autocorrelation function as a function of upper integrand limit:

\[y(t) = \int_0^{t}<P_{\alpha\beta}(0)P_{\alpha\beta}(\tau)> d\tau\]

In the ideal case this graph will converge to a constant value as t → \(\infty\). This usually requires very long simulations. Your graph may look different.

In this case, the integral seems to converge to a value of 1.03 × 10^-9 hartree² fs bohr^-6.

To calculate the viscosity, we also need the volume \(V = 13.333^3 = 2370\) ų.

Converting everything to SI units (where 1 hartree bohr^-3 = 2.942 × 10¹³ Pa), we get

\[\eta = \frac{2370 \times 10^{-30}}{1.38 \times 10^{-23} \times 298} \times 1.03 \times 10^{-9} \times 2.942 \times 10^{13} \times 2.942 \times 10^{13} \times 10^{-15}\]

\[\eta = 0.00051 \rm{\ Pa\ s} = 0.51 \rm{\ mPa\ s}\]

The experimental room-temperature viscosity of benzene is about 0.5 mPa s.

Important considerations

• The simulation needs to be well equilibrated

• The viscosity is often sensitive to the supercell size. Use a larger supercell with more benzene molecules for a better value.

• The integral must converge in the long time-limit. You may want to increase the maximum correlation time (above given as 3000 frames or 15 ps) to a larger value, and increase the simulation time.

• The biggest contribution to the integral comes at short correlation times. It is therefore necessary set the Sample frequency to a small number, ideally even smaller than the 5 fs we used here.

• In general, make sure to equilibrate the density of the liquid first using an NPT simulation, before running an NVT simulation as above. (the above example used the experimental density, but that may not be suitable for the used force field)

• Better and more reliable viscosity values can most likely be obtained using the APPLE&P or GFNFF force fields.

Python script

Both the MD simulation and the calculation of the autocorrelation function (ACF) integral can be scripted with our PLAMS scripting framework. However, extracting the converged viscosity value can be difficult to do automatically. We therefore recommend you to read off the value manually (and not simply use the output from the script).

• viscosity.py and the xyz files for 250 K, 300 K, 350 K.

The three xyz-files each contain 16 benzene molecules whose lattice constants have been equilibrated at 250, 300 and 350 K. The script can be run as is and will perform 3 MD runs followed by an analysis. It also provides a method to plot the viscosities so that you are able to read off the converged value yourself, however the script also tries to automatically guess the converged value. It is recommended however to read them off yourself.

Example output:

[18.02|13:41:55] PLAMS working folder: /home/hordijk/Viscosity/Benzene
Beginning viscosity calculations Benzene
I will perform 3 MD runs
Starting MD runs ...
[18.02|13:41:55] JOB Benzene_T250 STARTED
[18.02|13:41:55] JOB Benzene_T300 STARTED
[18.02|13:41:55] JOB Benzene_T350 STARTED
[18.02|13:41:55] Waiting for job Benzene_T250 to finish
[18.02|13:41:55] JOB Benzene_T250 RUNNING
[18.02|13:41:56] JOB Benzene_T300 RUNNING
[18.02|13:41:56] JOB Benzene_T350 RUNNING
[18.02|14:21:18] JOB Benzene_T300 FINISHED
[18.02|14:21:40] JOB Benzene_T300 SUCCESSFUL
[18.02|15:02:37] JOB Benzene_T250 FINISHED
[18.02|15:03:00] JOB Benzene_T250 SUCCESSFUL
[18.02|15:03:00] Waiting for job Benzene_T350 to finish
[18.02|15:42:53] JOB Benzene_T350 FINISHED
[18.02|15:43:04] JOB Benzene_T350 SUCCESSFUL
Starting analysis runs ...
[18.02|15:43:11] JOB Benzene_T250_analysis STARTED
[18.02|15:43:11] JOB Benzene_T300_analysis STARTED
[18.02|15:43:11] JOB Benzene_T350_analysis STARTED
[18.02|15:43:11] Waiting for job Benzene_T250_analysis to finish
[18.02|15:43:11] JOB Benzene_T250_analysis RUNNING
[18.02|15:43:11] JOB Benzene_T300_analysis RUNNING
[18.02|15:43:11] JOB Benzene_T350_analysis RUNNING
[18.02|15:43:24] JOB Benzene_T250_analysis FINISHED
[18.02|15:43:24] JOB Benzene_T350_analysis FINISHED
[18.02|15:43:24] JOB Benzene_T300_analysis FINISHED
[18.02|15:43:24] JOB Benzene_T250_analysis SUCCESSFUL
[18.02|15:43:24] JOB Benzene_T350_analysis SUCCESSFUL
[18.02|15:43:24] JOB Benzene_T300_analysis SUCCESSFUL
[18.02|15:43:25] PLAMS run finished. Goodbye
All calculations are now done
Temperature +- SD (K) | Viscosity (mPa s)
   249.59 +- 14.74  | 1.47419
   300.49 +- 18.50  | 0.47794
   349.92 +- 21.23  | 0.17027

Important

These integrals are not converged, and reading off the viscosities from these plots is unreliable! To get more reliable results, run longer simulations and set SamplingFreq to a smaller number.