Example 3: Setting up finite difference calculation for CO on Pd(111)

This example uses tools in coolvib/tools/ to perform finite difference displacements with the help of ASE for the tensor and normal mode case.

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 perform a finite-difference calculation using ASE and FHI-aims
to generate matrix elements for coolvib"""

import numpy as np
import os, sys
from ase.io import read
from ase.calculators.aims import Aims
import coolvib
from coolvib.tools.finite_difference import finite_difference
from coolvib.tools.mode_displacement import mode_displacement
from coolvib.routines.output import print_matrix

input1='start.in'
ab=read(input1)

#CALCULATOR
basisset = 'light'

#setting up Aims calculator in ASE
calc=Aims(xc='PBE',
    run_command = 'mpirun -np 4 <insert_path>/bin/aims.160210.mpi.x > aims.out',
    species_dir='<insert path>/aimsfiles/species_defaults/'+basisset,
    occupation_type = ['gaussian',0.1],
    sc_iter_limit = 100,
    #spin = 'collinear',
    relativistic = ['atomic_zora','scalar'],
    #default_initial_moment = 0,
    sc_accuracy_etot=1e-6,
    sc_accuracy_eev=0.001,
    sc_accuracy_rho=1e-5,
    sc_accuracy_forces=1e-3,
    load_balancing = True,
    #empty_states=10,
    k_grid = [12,12,1],
    restart_aims='wvfn.dat',
    output = ["eigenvectors",  "k_point_list",  "hamiltonian_matrix",  "overlap_matrix",  "band     0.0000    0.0000    0.0000    0.0000    0.0833    0.0000  2 ", 
    "band     0.0000    0.2500    0.0000    0.0000    0.3333    0.0000  2 ", 
    "band     0.0000    0.4167    0.0000    0.0000    0.5000    0.0000  2 ", 
    "band     0.0000    0.8333    0.0000    0.0833    0.0000    0.0000  2 ", 
    "band     0.0833    0.0833    0.0000    0.0833    0.1667    0.0000  2 ", 
    "band     0.0833    0.3333    0.0000    0.0833    0.4167    0.0000  2 ", 
    "band     0.0833    0.5000    0.0000    0.0833    0.6667    0.0000  2 ", 
    "band     0.0833    0.7500    0.0000    0.0833    0.8333    0.0000  2 ", 
    "band     0.1667    0.0000    0.0000    0.1667    0.0833    0.0000  2 ", 
    "band     0.1667    0.1667    0.0000    0.1667    0.3333    0.0000  2 ", 
    "band     0.1667    0.4167    0.0000    0.1667    0.5000    0.0000  2 ", 
    "band     0.1667    0.6667    0.0000    0.1667    0.7500    0.0000  2 ", 
    "band     0.1667    0.8333    0.0000    0.2500    0.0000    0.0000  2 ", 
    "band     0.2500    0.0833    0.0000    0.2500    0.1667    0.0000  2 ", 
    "band     0.2500    0.3333    0.0000    0.2500    0.4167    0.0000  2 ", 
    "band     0.2500    0.5000    0.0000    0.2500    0.6667    0.0000  2 ", 
    "band     0.2500    0.7500    0.0000    0.2500    0.8333    0.0000  2 ", 
    "band     0.3333    0.0000    0.0000    0.3333    0.0833    0.0000  2 ", 
    "band     0.3333    0.1667    0.0000    0.3333    0.3333    0.0000  2 ", 
    "band     0.3333    0.4167    0.0000    0.3333    0.5000    0.0000  2 ", 
    "band     0.3333    0.6667    0.0000    0.3333    0.7500    0.0000  2 ", 
    "band     0.3333    0.8333    0.0000    0.4167    0.0000    0.0000  2 ", 
    "band     0.4167    0.0833    0.0000    0.4167    0.1667    0.0000  2 ", 
    "band     0.4167    0.3333    0.0000    0.4167    0.4167    0.0000  2 ", 
    "band     0.4167    0.5000    0.0000    0.4167    0.6667    0.0000  2 ", 
    "band     0.4167    0.7500    0.0000    0.4167    0.8333    0.0000  2 ", 
    "band     0.5000    0.0000    0.0000    0.5000    0.0833    0.0000  2 ", 
    "band     0.5000    0.1667    0.0000    0.5000    0.3333    0.0000  2 ", 
    "band     0.5000    0.4167    0.0000    0.5000    0.5000    0.0000  2 ", 
    "band     0.5000    0.7500    0.0000    0.5833    0.0833    0.0000  2 ", 
    "band     0.5833    0.4167    0.0000    0.5833    0.7500    0.0000  2 ", 
    "band     0.6667    0.0833    0.0000    0.6667    0.4167    0.0000  2 ", 
    "band     0.6667    0.7500    0.0000    0.7500    0.0833    0.0000  2 ", 
    "band     0.7500    0.4167    0.0000    0.7500    0.7500    0.0000  2 ", 
    "band     0.8333    0.0833    0.0000    0.8333    0.4167    0.0000  2 ", 
    "band     0.8333    0.7500    0.0000    0.9167    0.0833    0.0000  2 ", 
    "band     0.9167    0.4167    0.0000    0.9167    0.7500    0.0000  2 ", 
    ], 
) 

ab.set_calculator(calc)

indices1=[4,5]
f = finite_difference(ab,indices=indices1, delta=0.0025)

f.run()
f.summary()
f.write_jmol()

#printing the Hessian
print 'Vibrational hessian just for reference'
print_matrix(f.modes)

fl = open('hessian', 'w')
for i in range(len(f.hnu)):
    string = ''
    modes = f.modes
    for j in range(len(f.hnu)):
        string += ' {0:14.8f} '.format(modes[i,j])
    fl.write(string+'\n')
fl.close()

######
# Now we have generated all input files necessary to run example1
#####

# assert 0

#Calculate finite difference along a given mode

print ' Now calculating finite difference along a given normal mode'
#follow the internal stretch mode

#You can define your own mode displacement or use one calculated from 
#a Hessian calculation. It is important to use displacements in 
#mass-weighted cartesian coordinates, 
#that means, e*sqrt(mass)

#this is the last mode ... internal stretch 
mode = np.array([
     [0.,0.,0.],
     [0.,0.,0.],
     [0.,0.,0.],
     [0.,0.,0.],
     [0.,0.,0.84925],
     [0.,0.,-0.528],
      ])
#We pick the right mode from the hessian
#in this case the internal stretch mode

#mode = f.get_mode(-1)
#we normalize
mode /= np.linalg.norm(mode)
print mode
#we mass-weight
# print np.diag(ab.get_masses())
# w = np.diag(np.sqrt(ab.get_masses()))
# print w
# mode = np.dot(mode,w)
# print mode

mode_displacement(ab, mode, name='mode', disp=0.0025)


Detailed Explanation

In the above code we first initialize an ASE atoms object and an ASE FHI-aims calculator object. The output string is essential to guarantee that all input data is being dumped to file for each displacement. The tool coolvib/tools/fhiaims_add_outputoptions_to_control.py helps to identify this complicated string.

After setting the calculator we initialize the finite_difference object and run it. The result is a Hessian calculation that also dumps the necessary friction input files into folders.

mode_displacement does the same for a given normal mode