Add STO example in the tests folder

tests/sprkkr/STO/STO.py

@@ -0,0 +1,146 @@
+# 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")