147 lines
5.3 KiB
Python
147 lines
5.3 KiB
Python
|
# coding: utf8
|
||
|
|
# vim: set et sw=4 ts=4 fdm=indent nu cc=+1:
|
||
|
|
|
||
|
|
from ase.lattice.tetragonal import SimpleTetragonalFactory
|
||
|
|
from ase.visualize import view
|
||
|
|
from ase.spacegroup import get_spacegroup
|
||
|
|
from sprkkr.calculator import SPRKKR
|
||
|
|
from msspec.calculator import MSSPEC, XRaySource
|
||
|
|
from msspec.utils import *
|
||
|
|
import numpy as np
|
||
|
|
import logging
|
||
|
|
import sys
|
||
|
|
import os
|
||
|
|
import glob
|
||
|
|
|
||
|
|
|
||
|
|
logging.basicConfig(level=logging.DEBUG)
|
||
|
|
logger = logging.getLogger(__name__)
|
||
|
|
|
||
|
|
|
||
|
|
# Define a Perovskite Factory class
|
||
|
|
class PerovskiteFactory(SimpleTetragonalFactory):
|
||
|
|
bravais_basis = [[0, 0, 0.0], [0.5, 0.5, 0.5],[0.5, 0, 0], [0, 0.5, 0],
|
||
|
|
[0.0, 0.0, 0.5]]
|
||
|
|
element_basis = (0, 1, 2, 2, 2)
|
||
|
|
|
||
|
|
ABO3 = Perovskite = PerovskiteFactory()
|
||
|
|
|
||
|
|
# Generate the base STO cell
|
||
|
|
a0 = 3.905
|
||
|
|
STO = Perovskite(('Sr', 'Ti', 'O'),
|
||
|
|
latticeconstant={'a': a0, 'c/a': 1.},
|
||
|
|
size=(1,1,1))
|
||
|
|
prefix = "SrTiO3"
|
||
|
|
|
||
|
|
if 'view' in sys.argv:
|
||
|
|
view(STO)
|
||
|
|
|
||
|
|
# ########## SPRKKR part
|
||
|
|
if 'sprkkr' in sys.argv:
|
||
|
|
# create a SPRKKR calculator
|
||
|
|
calc = SPRKKR(label=f"{prefix}/{prefix}")
|
||
|
|
|
||
|
|
# attach the atoms to the calculator object
|
||
|
|
calc.set_atoms(STO)
|
||
|
|
|
||
|
|
# Here is how to modify input file
|
||
|
|
# calc.input.control_section.set(DATASET="Fe", ADSI="SCF", POTFIL="Fe.pot")
|
||
|
|
# calc.input.tau_section.set(nktab=250)
|
||
|
|
# launch kkrscf
|
||
|
|
calc.get_potential_energy()
|
||
|
|
|
||
|
|
# calc2 = SPRKKR(label="Fe/Fe", task="phagen")
|
||
|
|
# calc2.set_atoms(Fe)
|
||
|
|
# calc2.input.control_section.set(POTFIL="Fe.pot_new")
|
||
|
|
# calc2.phagen()
|
||
|
|
|
||
|
|
#
|
||
|
|
# EXPORT POTENTIAL FOR PHAGEN
|
||
|
|
#
|
||
|
|
#change task and command
|
||
|
|
calc.set_command('PHAGEN')
|
||
|
|
# Change output file
|
||
|
|
calc.set_outfile(f"{prefix}_phagen.out")
|
||
|
|
#to change task we need to replace input file tempate
|
||
|
|
calc.set_inpfile(f"{prefix}_phagen.inp")
|
||
|
|
calc.input.control_section.set(DATASET="PHAGEN", ADSI="PHAGEN")
|
||
|
|
#set potetential file to converged potential
|
||
|
|
conv_potfile=os.path.join(calc.potfile+'_new')
|
||
|
|
|
||
|
|
calc.set_potfile(conv_potfile)
|
||
|
|
#run given task
|
||
|
|
calc.phagen()
|
||
|
|
|
||
|
|
# ######### MsSpec part
|
||
|
|
if 'msspec' in sys.argv:
|
||
|
|
###########################################################################
|
||
|
|
# Build the SrTiO3 cluster #
|
||
|
|
###########################################################################
|
||
|
|
# We start by loading the potential since it updates info in the STO cell
|
||
|
|
pot = SPRKKRPotential(STO, "SrTiO3/SrTiO3.pot_new",
|
||
|
|
*glob.glob("SrTiO3/*PHAGEN.pot"))
|
||
|
|
|
||
|
|
# tag each atom in the STO object to easily set the emitter in the
|
||
|
|
# hemispherical_cluster function
|
||
|
|
[atom.set('tag', ('Sr', 'Ti', 'O').index(atom.symbol)) for atom in STO]
|
||
|
|
|
||
|
|
# Create a hemispherical cluster centered on a Ti emitter from the
|
||
|
|
# STO primitive cell used in the SPRKKR step.
|
||
|
|
# For this example we use 4 planes and the Ti emitter (tag #1) is
|
||
|
|
# located in the 2nd plane (numbering starts at 0)
|
||
|
|
cluster = hemispherical_cluster(STO,
|
||
|
|
planes=4,
|
||
|
|
emitter_plane=1, emitter_tag=1)
|
||
|
|
|
||
|
|
# The created cluster is centered on the emitter atom, so defining
|
||
|
|
# the absorber attribute is straight forward:
|
||
|
|
cluster.absorber = get_atom_index(cluster, 0, 0, 0)
|
||
|
|
|
||
|
|
###########################################################################
|
||
|
|
# Set up the PhotoElectron Diffraction calculator #
|
||
|
|
###########################################################################
|
||
|
|
# Create the calculator
|
||
|
|
calc = MSSPEC(spectroscopy='PED', algorithm='expansion', folder='PED')
|
||
|
|
|
||
|
|
# minimalistic set of parameter for XPD scan
|
||
|
|
calc.set_atoms(cluster)
|
||
|
|
calc.source_parameters.energy = 700
|
||
|
|
calc.calculation_parameters.scattering_order = 2
|
||
|
|
|
||
|
|
# Run a polar scan with the default MuffinTin potential for Ti(2p3/2)
|
||
|
|
data = calc.get_theta_scan(level='2p3/2')
|
||
|
|
|
||
|
|
# Now we use the previously generated SPRKKR potential
|
||
|
|
# and run the same calculation with appending the data to the previous one
|
||
|
|
calc.tmatrix_parameters.potential = pot
|
||
|
|
data = calc.get_theta_scan(level='2p3/2', data=data)
|
||
|
|
|
||
|
|
# To better see the differences, plot the normalized signal on the
|
||
|
|
# same graph
|
||
|
|
# pick up previous data
|
||
|
|
theta, muffintin_cs, sprkkr_cs = (data[0].theta.copy(),
|
||
|
|
data[0].cross_section.copy(),
|
||
|
|
data[1].cross_section.copy())
|
||
|
|
# divide by the max
|
||
|
|
sprkkr_cs /= sprkkr_cs.max()
|
||
|
|
muffintin_cs /= muffintin_cs.max()
|
||
|
|
# add a new dataset with those values
|
||
|
|
dset = data.add_dset('comparison')
|
||
|
|
dset.add_columns(theta=theta, sprkkr=sprkkr_cs, muffintin=muffintin_cs)
|
||
|
|
# add a view for this dataset
|
||
|
|
view = dset.add_view('Comparison',
|
||
|
|
title=r'Comparison of XPD polar scans for Ti(2p3/2)',
|
||
|
|
xlabel=r'$\Theta (\degree$)',
|
||
|
|
ylabel='Normalized Signal (a.u.)')
|
||
|
|
view.select('theta', 'sprkkr', legend='With SPRKKR potential')
|
||
|
|
view.select('theta', 'muffintin', legend='With internal MuffinTin potential')
|
||
|
|
view.set_plot_options(autoscale=True)
|
||
|
|
|
||
|
|
# Pop up the final result
|
||
|
|
data.view()
|
||
|
|
|
||
|
|
|
||
|
|
if len(sys.argv) <= 1:
|
||
|
|
print("Please specify either 'sprkkr', 'msspec' keywords or both "
|
||
|
|
"of them on the command line")
|