Pietas aims to be a Platform for InElastic Tunneling At Surfaces. It is a first-principles post-processing code based on VASP electronic structure calculation. We used finite difference method to calculate the change in wavefunctions along a vibration normal mode, according to the poineering works by Prof. N. Lorente et al. [PRL, 85, 2997(2000); Farady Discuss. 117, 277(2000)]. The current implementation is baed on our previous FORTRAN code completed around 2008 [H. Ren, et al. JCP, 130, 134707(2009)]. Currently only the Tersorff-Hamann approximation is implemented.

• A modern operating system that can run CPython with the following modules...
• Python (2.7+ or 3.4+)
• numpy (1.9.0+)
• scipy (0.15.0+)
• VASP (4.6+ or 5.3+)

Unpack the package to any directory you can access (say $PIETAS), make sure all the required modules are available in the python environment you are using. Optionaly you can also create a shortcut for your lcoal environment:

$ cd pietas
$ export PIETAS=`pwd`

I personaly recommend the Anacoda Python distribution, which contains most of the modules widely used in the scientific computation community, as well as Intel MKL acceleration with auto parallization via openmp.

The examples lies in the directory $PIETAS/examples. Currently we have only one classic case CO adsorbed on Cu(100) surface ($PIETAS/examples/CO\@Cu100).

All the input files reqired for the VASP electronic structure calculation and frequency analysis are available in this directory.

At the very beginning, we will perform a VASP frequency calculation with IBRION = 5. A step size of POTIM = 0.02 Angstrom is used as suggested in the VASP 4.6 manual (the default value has been changed to 0.015 if unreasonable large values set in INCAR). We use a 4x4x1 k-sampling in the 1st Brillouin zone centerd at the Gamma point (the KPOINT file). PBE functional is used to describe the exchange-correlation effects, and the core electrons are described by the PAW scheme with frozen core approximation. There's also a POSCAR file prepared for the calculation, which is optimized with the same paramerter by fixing the bottom two Cu layers. An OUTCAR.comp.tgz file is also attached for benchmark purpose.

1. Frequency analysis

   Now you can perform the frequency calculation with your VASP executable, or you can also just unpack OUTCAR.comp.tgz for a usable OUTCAR, which is required by the genmode module disscussed in the next step. We also provide a tool for visualize the vibrational modes by using the program MOLDEN. The following command

   $ python viewmode.py -i /path/to/freq/OUTCAR

   will generate the vibrational modes in the MOLDEN format, then you can select a vibrational mode by its frequency (in meV) and visualize it. Please note that the OpenGL version (gmolden) of MOLDEN might perform better on modern hardwares.

2. Generate dispalced geometry configurations

   This step will generate the displaced geometry configurations along the vibrational modes specified. Since the displacement dQ is respect to the mean squared amplitude for each mode, which would have a very small value for high frequency vibrations, we use a scale factor to enlongate/compress dQ to ensure the real displacements are large enough beyond the error bar of DFT numerical calculation. By default this scale factor is selected as 1, which is somewhat reasonable for vibration modes with moderate frequency, such as the two-fold degenerated frustrated rotation mode of CO@Cu(100) at ~31 meV.

   To generate the displaced structures, just use the command:

   $ python genmod.py -i /path/to/freq/OUTCAR

   the displaced structures will be created and saved in the directory ./f1.0, with file names POSCAR with suffixes .+/-mddd, where the plus/minus sign denotes forward/backward displacement. The character m is followed by a three-digit number corresponds to the ddd-th vibration mode. Displaced structures correspond to all the vibration modes will be generated by default.

   You can also specify only one mode to be generated:

   $ phthon genmode.py -i /path/to/freq/OUTCAR -m 3

   The scaling factor can be specified by the argument -f:

   $ python genmode.py -i /path/to/freq/OUTCAR -f 2.0

3. Static electronic structure calculation

   Once the displaced structures are ready, we can perform static calculations for them. We need two static calculations for each mode, one for backward displacement, and the other for forward. In the directory ./examples/CO\@Cu100/f1.0/mode3/, we prepared two subdirectories ./b and ./f for these two opposite displacements. Just copy the generated POSCAR.+/-mddd file to the corresponding directory and rename to POSCAR. The PAW dataset POTCAR is exactly the same with that used in step 1, so just make a copy.

   Make sure that you are preparing a static calculation with IBRION=-1, then make a VASP run.

4. IETS calculation

   The IETS calculation uses the wavefunctions stored in the WAVECAR files generated by the static calculations in the last step. For each mode, we need three WAVECAR files, one for the equilibrium configuration, and the others for backward/forward displaced structures.

   The pietas module is controled by a one-file input, with the format YAML, which is self-explained by comments. You can find an input file named f1.0-mode3.yml in the $PIETAS/examples/CO\@Cu100/f1.0 directory. To run the pietas calculation for mode3, a slight modification of this input files might be necessary.

1. get into the working directory
    
    ```shell
    cd $PIETAS
    ```
2. copy the input file into the working directory

    ```shell
    cp examples/CO\@Cu100/f1.0/f1.0-mode3.yml .
    ```
3. edit the input file to suit your enviroment, mostly the paths of `WAVECAR` and `OUTCAR` files
    - the `wavecar` block:

        modify the paths of the three `WAVECAR` files as in your system, for example:
        ```yaml
        backward    :    ./examples/CO\@Cu100/f1.0/mode3/b/WAVECAR
        ```
    + the `outcar` block:

        modify the path to point to a complete `OUTCAR` file generated with the equilibrium structure.
    + the `output path` block

        modify to pointing to a suitable path you want to store the results.

Then run the pietas calculation:

$ python pietas.py -i f1.0-mode3.yml

The calculation will take several minutes depending on your hardware and parameters (such as sigma, cutoff, and real grid density)

Pietas will loop over spin and k-points, the change in LDOS (principal and inelastic) and the LDOS itself will be stored in the output path. The total or k-point weighted LDOS and its change are also saved in the same directory. All the output files are now stored in the npy format which is native to numpy, and are ready to use in matplotlib.

• Complete output generation, at least for XCrysDen (XSF) and NetCDF (nc) formats.
• Add all-electron support for the PAW method.
• Add Bardeen approximation (with tip).
• Add WKB approximation.