Tribology

Tribology is the study of two interfaces sliding relative to one another, as encountered in frictional sliding or adhesion. Molecular dynamics simulations of representative volume elements of tribological interfaces are routinely used to gain insights into the atomistic mechanisms underlying friction and wear. The module matscipy.pressurecoupling provides tools to perform such simulations under a constant normal load and sliding velocity.

The example below shows how to perform an initial fast pressure equilibration of an interface prior to sliding. Afterwards, during sliding, we apply the pressure coupling by Pastewka et al. to dynamically adjust the distance between the two surfaces according to the local pressure. The algorithm ensures mechanical boundary conditions that account for the inertia of the bulk material which is not explicitly included in the simulation.

System setup

Let’s first create an exemplary sliding interface consisting of two silicon crystals:

from ase.lattice.cubic import Diamond

# create two slabs
slab1 = Diamond(directions=[[1, 0, 0], [0, 1, 0], [0, 0, 1]], symbol='Si',
                pbc=(1, 1, 1), latticeconstant=5.431, size=(2, 2, 2))

slab2 = slab1.copy()

# merge them
slab2.translate([1, 1, slab1.get_cell()[2][2] + 1.2])
atoms = slab1 + slab2

# remove pbc along z direction
atoms.center(axis=2, vacuum=0)
atoms.set_pbc([True, True, False])

Initial pressure equilibration

To apply a normal load along the z direction to this system, we will fix the two lowest atomic layers of the lower crystal. The corresponding atoms are defined with the mask (bottom_mask). The two topmost atomic layers of the upper crystal will be treated rigidly (top_mask). The z position of this rigid region can adapt according to the normal pressure. To reuse the masks in subsequent simulation runs (sliding or restarts), we save them in files.

import numpy as np

bottom_mask = atoms.get_positions()[:, 2] < 1.4  # mask for lower fixed atoms
top_mask = atoms.get_positions()[:, 2] > atoms.get_cell()[2][2] - 1.4  # mask for upper rigid atoms

np.savetxt("bottom_mask.txt", bottom_mask)
np.savetxt("top_mask.txt", top_mask)

We now specify the numerical parameters of the MD simulation. A Langevin thermostat will be used in the y direction along which neither the normal load nor the sliding motion are applied. For simplicity, we will thermalize all atoms which are not part of the constraint regions. This makes sense for small systems which cannot have a dedicated thermostat region.

from ase.units import GPa, fs

dt = 5.0 * fs  # MD time step
C11 = 166.0 * GPa  # material constant
M_factor = 1.0  # scaling factor for lid mass during equilibration,
# 1.0 will give fast equilibration for expensive calculators

Pdir = 2  # index of cell axis along which the normal pressure is applied
P = 1.0 * GPa  # target normal pressure
v = 0.0  # no sliding yet, only apply pressure
vdir = 0  # index of cell axis along which sliding happens
T = 300.0  # target temperature for thermostat in K

t_langevin = 75.0 * fs  # time constant for Langevin thermostat
gamma_langevin = 1. / t_langevin  # derived Langevin parameter
t_integrate = 100.0 * fs  # simulation time
steps_integrate = int(t_integrate / dt)  # number of simulation steps

Next, we set up the calculation by attaching the constraint and a calculator to the atoms object. For the initial pressure equilibration, we critically damp the motion of the upper rigid layer while not increasing its mass. This is useful for a fast pressure equilibration. Interatomic interactions are modelled with the Stillinger-Weber potential for silicon. For a fast temperature convergence, since we typically start from a local minimum, we also set the initial temperature of the system to twice the target temperature.

from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber, Stillinger_Weber_PRB_31_5262_Si
from matscipy.calculators.manybody import Manybody
from matscipy import pressurecoupling as pc
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution

damp = pc.FixedMassCriticalDamping(C11, M_factor)  # specifies the critical damping with a fixed mass 
slider = pc.SlideWithNormalPressureCuboidCell(top_mask, bottom_mask, Pdir,
                                              P, vdir, v, damp)  # ASE constraint for pressure coupling

atoms.set_constraint(slider)  # attach constraint to atoms object

MaxwellBoltzmannDistribution(atoms, temperature_K=2 * T)  # initialize temperature

# clear momenta in constraint regions, otherwise upper rigid region might run away
atoms.arrays['momenta'][top_mask, :] = 0
atoms.arrays['momenta'][bottom_mask, :] = 0

calc = Manybody(**StillingerWeber(Stillinger_Weber_PRB_31_5262_Si))  # specify calculator
atoms.set_calculator(calc)  # attach calculator

/usr/lib/python3/dist-packages/scipy/__init__.py:146: UserWarning: A NumPy version >=1.17.3 and <1.25.0 is required for this version of SciPy (detected version 1.26.4
  warnings.warn(f"A NumPy version >={np_minversion} and <{np_maxversion}"

/tmp/ipykernel_9020/829096581.py:19: DeprecationWarning: Please use atoms.calc = calc
  atoms.set_calculator(calc)  # attach calculator

Finally, we setup the integrator and run the simulation. This will create an ASE trajectory file of the MD run and a log file that tracks the current status of the simulation in terms of the system temperature, the system height, the shear stress etc. The log file can be read using the class (SlideLog).

from ase.io import Trajectory
from ase.md.langevin import Langevin
from io import open

# only thermalize middle region in one direction
temps = np.zeros((len(atoms), 3))
temps[slider.middle_mask, slider.Tdir] = T
gammas = np.zeros((len(atoms), 3))
gammas[slider.middle_mask, slider.Tdir] = gamma_langevin

# set up integrator
integrator = Langevin(atoms, dt, temperature_K=temps,
                      friction=gammas, fixcm=False)

# set up trajectory file
trajectory = Trajectory('equilibrate_pressure_01.traj', 'w', atoms) 

# set up logger
log_handle = open('log_equilibrate_01.txt', 'w', 1, encoding='utf-8')  # 1 means line buffered,
logger = pc.SlideLogger(log_handle, atoms, slider, integrator)

logger.write_header()  # write header of log file
integrator.attach(logger)  # attach logger to integrator
integrator.attach(trajectory)  # attach trajectory to integrator
integrator.run(steps_integrate)  # run the simulation
log_handle.close()
trajectory.close()

Restarting a pressure equilibration

To restart the pressure equilibration, we start by specifying the parameters of the MD simulation, read the previous trajectory to initialize the current status of the system, and read the masks for the constraint regions.

from ase.units import GPa, fs
from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber, Stillinger_Weber_PRB_31_5262_Si
from matscipy.calculators.manybody import Manybody
from matscipy import pressurecoupling as pc
from ase.io import read
import numpy as np

dt = 5.0 * fs  # MD time step
C11 = 166.0 * GPa  # material constant
M_factor = 1.0  # scaling factor for lid mass during equilibration

Pdir = 2  # index of cell axis along normal pressure is applied
P = 1.0 * GPa  # target normal pressure
v = 0.0  # no sliding yet, only apply pressure
vdir = 0  # index of cell axis along sliding happens
T = 300.0  # target temperature for the thermostat

t_langevin = 75.0 * fs  # time constant for Langevin thermostat
gamma_langevin = 1. / t_langevin  # derived Langevin parameter
t_integrate = 100.0 * fs  # simulation time
steps_integrate = int(t_integrate / dt)  # number of simulation steps

# get atoms from trajectory to also initialize correct velocities
atoms = read('equilibrate_pressure_01.traj')

bottom_mask = np.loadtxt("bottom_mask.txt").astype(bool)
top_mask = np.loadtxt("top_mask.txt").astype(bool)

damp = pc.FixedMassCriticalDamping(C11, M_factor)
slider = pc.SlideWithNormalPressureCuboidCell(top_mask, bottom_mask,
                                              Pdir, P, vdir, v, damp)
atoms.set_constraint(slider)

calc = Manybody(**StillingerWeber(Stillinger_Weber_PRB_31_5262_Si))  # specify calculator
atoms.set_calculator(calc)

/tmp/ipykernel_9020/3486272778.py:35: DeprecationWarning: Please use atoms.calc = calc
  atoms.set_calculator(calc)

Next, we again initialize the integrator.

from ase.md.langevin import Langevin

temps = np.zeros((len(atoms), 3))
temps[slider.middle_mask, slider.Tdir] = T
gammas = np.zeros((len(atoms), 3))
gammas[slider.middle_mask, slider.Tdir] = gamma_langevin
integrator = Langevin(atoms, dt, temperature_K=temps,
                      friction=gammas, fixcm=False)

We specify the trajectory file and the log file and run the simulation. Note that ASE automatically writes the initial step of a run, i.e., if you choose to append to the previous trajectory and log files, both will contain the information about the initial configuration twice, since this configuration has already been written to the files at the end of the previous run.

from ase.io import Trajectory
from io import open

trajectory = Trajectory('equilibrate_pressure_02.traj', 'w', atoms)

with open('log_equilibrate_01.txt', 'r', encoding='utf-8') as log_handle:
    step_offset = pc.SlideLog(log_handle).step[-1]  # read last step

log_handle = open('log_equilibrate_02.txt',
                  'w', 1, encoding='utf-8')  # line buffered
logger = pc.SlideLogger(log_handle, atoms, slider, integrator, step_offset)
logger.write_header()  # write header of log file

integrator.attach(logger)
integrator.attach(trajectory)
integrator.run(steps_integrate)
log_handle.close()
trajectory.close()

Sliding simulation

After the pressure equilibration, we can start to apply a shear deformation to the system by applying a constant sliding velocity to the upper rigid layer. First, we specify the parameters of the MD simulation, read the equilibration trajectory to initialize the current status of the system, and read the masks for the constraint regions.

from ase.units import GPa, fs, m, s
from ase.io import read
import numpy as np
from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber, Stillinger_Weber_PRB_31_5262_Si
from matscipy.calculators.manybody import Manybody

dt = 5.0 * fs  # MD time step
C11 = 166.0 * GPa  # material constant
p_c = 0.20  # empirical cutoff parameter value to remove high-frequency oscillations of the system 

Pdir = 2  # index of cell axis along normal pressure is applied
P = 1.0 * GPa  # target normal pressure
v = 100.0 * m / s  # constant sliding speed
vdir = 0  # index of cell axis along sliding happens
T = 300.0  # target temperature for thermostat

t_langevin = 75.0 * fs  # time constant for Langevin thermostat
gamma_langevin = 1. / t_langevin  # derived Langevin parameter
t_integrate = 100.0 * fs  # simulation time
steps_integrate = int(t_integrate / dt)  # number of simulation steps

# get atoms from trajectory to also initialize correct velocities
atoms = read('equilibrate_pressure_02.traj')

bottom_mask = np.loadtxt("bottom_mask.txt").astype(bool)
top_mask = np.loadtxt("top_mask.txt").astype(bool)

calc = Manybody(**StillingerWeber(Stillinger_Weber_PRB_31_5262_Si))  # specify calculator
atoms.set_calculator(calc)

/tmp/ipykernel_9020/1365514820.py:29: DeprecationWarning: Please use atoms.calc = calc
  atoms.set_calculator(calc)

Since we change the pressure coupling algorithm, we need to reset the velocities of the upper rigid layer.

velocities = atoms.get_velocities()
velocities[top_mask, Pdir] = 0.0

atoms.set_velocities(velocities)

In analogy to the pressure equilibration, we set up the constraint for the simulation, but this time we use the pressure coupling algorithm of Pastewka et al. rather than critical damping with a fixed mass.

from matscipy import pressurecoupling as pc

damp = pc.AutoDamping(C11, p_c)  # Damping by Pastewka et al.
slider = pc.SlideWithNormalPressureCuboidCell(top_mask, bottom_mask,
                                              Pdir, P, vdir, v, damp)
atoms.set_constraint(slider)

Afterwards, in analogy to the pressure equilibration, we initialize the integrator, the trajectory and log files, and run the simulation.

from ase.io import Trajectory
from ase.md.langevin import Langevin
from io import open

temps = np.zeros((len(atoms), 3))
temps[slider.middle_mask, slider.Tdir] = T
gammas = np.zeros((len(atoms), 3))
gammas[slider.middle_mask, slider.Tdir] = gamma_langevin
integrator = Langevin(atoms, dt, temperature_K=temps,
                      friction=gammas, fixcm=False)

trajectory = Trajectory('slide_01.traj', 'w', atoms)
log_handle = open('log_slide_01.txt', 'w', 1, encoding='utf-8')  # line buffered
logger = pc.SlideLogger(log_handle, atoms, slider, integrator)
# log can be read using pc.SlideLog (see docstring there)
logger.write_header()

integrator.attach(logger)
integrator.attach(trajectory)
integrator.run(steps_integrate)
log_handle.close()
trajectory.close()

Restarting a sliding simulation

To restart the sliding simulation, the following script could be used.

from ase.io import Trajectory, read
from ase.units import GPa, fs, m, s
import numpy as np
from ase.md.langevin import Langevin
from matscipy import pressurecoupling as pc
from io import open
from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber, Stillinger_Weber_PRB_31_5262_Si
from matscipy.calculators.manybody import Manybody

# Parameters
dt = 5.0 * fs  # MD time step
C11 = 166.0 * GPa  # material constant
p_c = 0.20  # empirical cutoff parameter value to remove high-frequency oscillations of the system

Pdir = 2  # index of cell axis along which normal pressure is applied
P = 1.0 * GPa  # target normal pressure
v = 100.0 * m / s  # constant sliding speed
vdir = 0  # index of cell axis along which sliding happens
T = 300.0  # target temperature for thermostat

t_langevin = 75.0 * fs  # time constant for Langevin thermostat
gamma_langevin = 1. / t_langevin  # derived Langevin parameter
t_integrate = 100.0 * fs  # simulation time
steps_integrate = int(t_integrate / dt)  # number of simulation steps

# get atoms from trajectory to also initialize correct velocities
atoms = read('slide_01.traj')

bottom_mask = np.loadtxt("bottom_mask.txt").astype(bool)
top_mask = np.loadtxt("top_mask.txt").astype(bool)

# set up sliding constraints
damp = pc.AutoDamping(C11, p_c)
slider = pc.SlideWithNormalPressureCuboidCell(top_mask, bottom_mask,
                                              Pdir, P, vdir, v, damp)
atoms.set_constraint(slider)

# set up calculator
calc = Manybody(**StillingerWeber(Stillinger_Weber_PRB_31_5262_Si))
atoms.set_calculator(calc)

# set up integrator
temps = np.zeros((len(atoms), 3))
temps[slider.middle_mask, slider.Tdir] = T
gammas = np.zeros((len(atoms), 3))
gammas[slider.middle_mask, slider.Tdir] = gamma_langevin
integrator = Langevin(atoms, dt, temperature_K=temps,
                      friction=gammas, fixcm=False)

# specify log file and trajectory file and run the simulation
trajectory = Trajectory('slide_02.traj', 'w', atoms)

with open('log_slide_01.txt', 'r', encoding='utf-8') as log_handle:
    step_offset = pc.SlideLog(log_handle).step[-1]  # read last step

log_handle = open('log_slide_02.txt', 'w',
                  1, encoding='utf-8')  # line buffered
logger = pc.SlideLogger(log_handle, atoms, slider, integrator, step_offset)
logger.write_header()

integrator.attach(logger)
integrator.attach(trajectory)
integrator.run(steps_integrate)
log_handle.close()
trajectory.close()

/tmp/ipykernel_9020/2714859955.py:40: DeprecationWarning: Please use atoms.calc = calc
  atoms.set_calculator(calc)

To remove all simulation files created during this tutorial, you can use the following code.

import os

for tmp_file in ['log_equilibrate_01.txt', 'log_equilibrate_02.txt',
                 'log_slide_01.txt', 'log_slide_02.txt',
                 'equilibrate_pressure_01.traj', 'equilibrate_pressure_02.traj',
                 'slide_01.traj', 'slide_02.traj',
                 'top_mask.txt', 'bottom_mask.txt']:
    try:
        os.remove(tmp_file)
    except FileNotFoundError:
        continue