Example 1: Calculating a friction tensor from FHI-aims data: CO on Cu(100)¶
In this example we calculate the friction tensor from data previously calculated with FHI-Aims. We do this for CO adsorbed at a top site of a Pd(111) surface. The corresponding friction tensor is a (6x6) matrix.
In order to run this example you need to download the FHI-Aims dataset from following URL and extract it in the example folder
Code¶
# This file is part of coolvib
#
# coolvib is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# coolvib is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with coolvib. If not, see <http://www.gnu.org/licenses/>.
"""
We calculate the (6x6) friction tensor for
a CO molecule adsorbed on a Pd(111) top site.
"""
import numpy as np
from ase.io import read
import os, sys
import coolvib
from scipy import linalg as LA
#from ase.visualize import view
#####DEFINE SYSTEM and PARAMETERS######
system = read('CO_on_Pd111/eq/geometry.in')
cell = system
active_atoms = [4,5] # meaning we have two atoms - C - 0 and O - 1
#view(system)
model = coolvib.workflow_tensor(system, code='aims', active_atoms=active_atoms)
finite_difference_incr = 0.0025
keywords = {
'discretization_type' : 'gaussian',
'discretization_broadening' : 0.01,
'discretization_length' : 0.01,
'max_energy' : 3.00,
'temperature' : 300,
'delta_function_type': 'gaussian',
'delta_function_width': 0.60,
'perturbing_energy' : 0.0,
'debug': 0,
}
print 'workflow initialized and keywords set'
######READ QM INPUT DATA###
model.read_input_data(
spin=False,
path='./CO_on_Pd111/',
filename='aims.out',
active_atoms=active_atoms,
incr=finite_difference_incr,
debug=0,
)
print 'successfully read QM input data'
############
#We can either calculate the spectral function and evaluate friction
#from there (discretization commands)
######CALCULATE SPECTRAL FUNCTION###
# model.calculate_spectral_function(mode='default', **keywords)
# # model.read_spectral_function()
# print 'successfully calculated spectral_function'
# model.print_spectral_function('nacs-spectrum.out')
# model.calculate_friction_tensor_from_spectrum(**keywords)
# model.plot_spectral_function()
#######CALCULATE FRICTION TENSOR###
#or we directly calculate friction
model.calculate_friction_tensor( **keywords)
print 'successfully calculated friction tensor'
model.analyse_friction_tensor()
model.print_jmol_friction_eigenvectors()
######project lifetimes along vibrations
hessian_modes = np.loadtxt('CO_on_Pd111/hessian')
for i in range(len(hessian_modes)):
hessian_modes[i,:]/=np.linalg.norm(hessian_modes[i,:])
from coolvib.constants import time_to_ps
ft= np.dot(hessian_modes,np.dot(model.friction_tensor,hessian_modes.transpose()))/time_to_ps
E, V = np.linalg.eigh(ft)
coolvib.routines.output.print_matrix(ft)
print 'Relaxation rates along vibrational modes in 1/ps'
for i in range(len(ft)):
print ft[i,i]
print 'Vibrational lifetimes along vibrational modes in 1/ps'
for i in range(len(ft)):
print 1./ft[i,i]
Detailed Explanation¶
System Setup¶
We start by setting up the system using an ASE atoms object and we initialize our workflow_tensor class with
import coolvib
model = coolvib.workflow_tensor(system, code='aims', active_atoms=[4,5])
, where active_atoms specifies the indices of the atoms we want to include in
the calculation of the friction tensor. Model is an instance of the coolvib.convenience.workflow_tensor.workflow_tensor
class that guides us through the calculation and bundles all data and functions.
Reading input data¶
After defining all keywords we can read the FHI-Aims data by
using coolvib.convenience.workflow_tensor.workflow_tensor.read_input_data()
.
The parser assumes that the input data is given as a set of directories with eq refering to the calculation at equilibrium position and a3c0+ or a3c0- refering to displacements in positive and negative directions for atom3 and the cartesian x component (c0 = x, c1 = y, c2 = z).
Calculating the spectral function¶
This is an optional step when calculating the friction tensor and at the moment only serves analysis purposes.
Using coolvib.convenience.workflow_tensor.workflow_tensor.calculate_spectral_function()
we can calculate the spectral functions along all cartesian components and all combinations of cartesian components. This includes loops over all possible excitations and construction of the electron-phonon spectral function on a discrete grid. This discretization is determined by following keywords:
- discretization_type
- Specifies which broadening function is used to discretize the individual excitations, for more details see
coolvib.routines.discretize_peak()
- discretization_broadening
- Specifies the width that is used to generate a finite width discrete representation of the peak and intensity (in eV)
- discretization_length
- Specifies the bin width of equidistant grid on which the spectral function is represented (in eV).
- max_energy
- Specifies the maximum extend beyond the fermi level for which the spectral function is calculated. The default is 3 eV. A safe value is 5 times the delta_function_width+perturbing_energy
The calculated spectral functions can be printed to file using coolvib.convenience.workflow_tensor.workflow_tensor.print_spectral_function()
. This file can be post-processed and further analyzed using scripts from coolvib.tools
.
Calculating the friction tensor and analysis¶
In order to finally arrive at the friction tensor one either has to integrate the spectral functions
over a delta function (see coolvib.routines.evaluate_delta_function()
). This is where the
main approximations of the underlying method are applied. The width, position and shape of the
delta function are controlled by following keywords:
- delta_function_type
- Specifies the type of delta function representation to be used
- delta_function_width
- Specifies the width of the delta function representation (in eV)
- perturbing_energy
- Specifies the position of the delta function representation. Typically this is chosen to be at the Fermi level under the quasi-static approximation (see Theory)
Convergence with respect to the decay rates and lifetimes has to be validated for the above parameters.
The function coolvib.convenience.workflow_tensor.workflow_tensor.calculate_friction_tensor_from_spectrum()
performs this integration and sets up the tensor.
One can directly calculate the friction tensor by using coolvib.convenience.workflow_tensor.calculate_friction_tensor()
, as it is done in this example
Analysis and printing of the tensor, its eigenvalues and eigenfunctions as well as the inverse of the eigenvalues (the principal lifetimes) can be done using coolvib.convenience.workflow_tensor.workflow_tensor.analyse_friction_tensor()
.
Results¶
The expected results are
Friction Tensor
0 1 2 3 4 50 3.2947e-01 -5.1951e-05 -1.4243e-04 -9.7699e-02 2.4844e-05 -1.1275e-04 1 -5.1951e-05 3.2955e-01 -1.4646e-04 1.2838e-05 -9.7716e-02 -1.1241e-04 2 -1.4243e-04 -1.4646e-04 1.6286e-01 4.2497e-05 4.2755e-05 -1.1424e-02 3 -9.7699e-02 1.2838e-05 4.2497e-05 2.9649e-02 -5.7152e-06 2.9358e-05 4 2.4844e-05 -9.7716e-02 4.2755e-05 -5.7152e-06 2.9654e-02 3.0400e-05 5 -1.1275e-04 -1.1241e-04 -1.1424e-02 2.9358e-05 3.0400e-05 6.7262e-03
- Friction Eigenvalues in 1/ps
- 0 1 2 3 4 5
0.0006 0.0006 0.0059 0.1637 0.3585 0.3586
Friction Eigenvectors
0 1 2 3 4 50 -2.7421e-01 -7.6918e-02 -1.3190e-04 6.8885e-04 -8.5161e-01 4.4006e-01 1 7.6930e-02 -2.7419e-01 -1.8838e-04 7.0623e-04 -4.4007e-01 -8.5161e-01 2 -3.7931e-05 -5.5541e-05 -7.2589e-02 9.9736e-01 1.0090e-03 3.3623e-04 3 -9.2294e-01 -2.5899e-01 8.2028e-04 -2.0170e-04 2.5301e-01 -1.3073e-01 4 2.5899e-01 -9.2295e-01 6.3231e-04 -2.1305e-04 1.3071e-01 2.5300e-01 5 -5.7038e-04 -7.3214e-04 -9.9736e-01 -7.2590e-02 4.1327e-04 1.3106e-04
- Principal lifetimes in ps
- 0 1 2 3 4 5
1605.4274 1597.3460 169.6452 6.1092 2.7896 2.7885