#!/usr/bin/env python
# encoding: utf-8
##############################################################################
# Author: Liam Deacon #
# #
# Contact: liam.deacon@diamond.ac.uk #
# #
# Copyright: Copyright (C) 2013-2014 Liam Deacon #
# #
# License: MIT License #
# #
# Permission is hereby granted, free of charge, to any person obtaining a #
# copy of this software and associated documentation files (the "Software"), #
# to deal in the Software without restriction, including without limitation #
# the rights to use, copy, modify, merge, publish, distribute, sublicense, #
# and/or sell copies of the Software, and to permit persons to whom the #
# Software is furnished to do so, subject to the following conditions: #
# #
# The above copyright notice and this permission notice shall be included in #
# all copies or substantial portions of the Software. #
# #
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR #
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, #
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL #
# THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER #
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING #
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER #
# DEALINGS IN THE SOFTWARE. #
# #
##############################################################################
"""
**model.py**
Provides convenience functions for generating input and calculating
atomic charge densities for use with the Barbieri/Van Hove phase
shift calculation package.
"""
from __future__ import print_function
from __future__ import division
import errno
import os
from copy import deepcopy
from shutil import move
from glob import glob
import phaseshifts.atorb as atorb
import phaseshifts.elements as elements
from phaseshifts.lib import libphsh
[docs]
class Atom(object):
"""
Atom class for input into cluster model for muffin-tin potential
calculations.
"""
def __init__(self, element, coordinates=[0.0, 0.0, 0.0], **kwargs):
"""
Constructor for Atom class.
Parameters
----------
element : str or int
This is either the elemental symbol, name or atomic number.
coordinates : list[float, float, float] or ndarray
The fractional coordinates within the unitcell in terms of the
basis vector a.
tag : str, optional
Add a name tag to this element (useful if dealing with multiple
atoms of the same element in a given model). Default is the
symbol for that element - numeric ids may be appended in the model
class.
radius : float, optional
The muffin-tin radius of the atom in Angstroms (default is to
lookup 'atmrad' in the element dictionary).
valence : int, optional
The valency of the atom (default is to assume neutral atom).
occupancy : float, optional
The fractional occupancy of the atom in the given site.
"""
self.element = elements.ELEMENTS[element]
self.coordinates = None # dummy
self.set_coordinates(coordinates)
self.name = self.element.name.title()
self.tag = self.element.symbol.title()
self.Z = self.element.protons
self.radius = self.element.atmrad
self.valence = 0
if "valence" in kwargs:
if "radius" not in kwargs:
# assume covrad for non-zero valency
self.radius = self.element.covrad
self.__dict__.update(kwargs)
# checks whether two atoms are equal w.r.t. name, radius and valence
def __eq__(self, other):
if isinstance(other, Atom):
return self.name == other.name and self.radius == other.radius and self.valence == other.valence
else:
return False
# checks whether two atoms are not equal w.r.t. name, radius and valence
def __neq__(self, other):
return not self.__eq__(other)
# reprinting of Atom object
def __repr__(self):
return str("Atom(%s, tag='%s', radius=%s, " "valence=%s)") % (
self.name,
self.tag,
self.radius,
self.valence,
)
# redefine hash method for checking uniqueness of class instance
def __hash__(self):
return hash(self.__repr__())
# set coordinates of atom within unitcell in terms of a
[docs]
def set_coordinates(self, coordinates):
self.coordinates = coordinates
self._coordinates = [r / 0.529 for r in coordinates]
# set valence of atom
[docs]
def set_valence(self, valency):
"""Sets the valency of the atom"""
self.valence = float(valency)
# set muffin-tin radius of atom
[docs]
def set_mufftin_radius(self, radius):
"""
Sets the muffin-tin radius of the atom in Angstroms.
"""
try:
self.radius = float(radius)
self._radius = self.radius / 0.529 # in Bohr radii
except (TypeError, ValueError):
return
[docs]
class Unitcell(object):
"""
Unitcell class
"""
def __init__(self, a, c, matrix_3x3, **kwargs):
"""
Constructor for the Unitcell class
Parameters
----------
a : float
The in-plane lattice vector in Angstroms
c : float
The out-of-plane lattice vector in Angstroms. For cubic systems this
will be equal to a.
matrix_3x3: ndarray
A 3x3 matrix describing the x,y,z construction of a,b,c basis vectors
of the unitcell. Units for x, y & z should be in terms of fractional
coordinates.
alpha : float, optional
Angle alpha in degrees.
beta : float, optional
Angle beta in degrees.
gamma : float, optional
Angle gamma in degrees.
"""
# Convert Angstrom input to Bohr radii
self.set_a(a)
self.set_c(c)
# Set basis vectors
self.set_vectors(matrix_3x3)
# Set xstal system
self.alpha = 90.0
self.beta = 90.0
self.gamma = 90.0
# Update additional information
self.__dict__.update(kwargs)
# checks if two class instances are equal w.r.t. name, radius & valence
def __eq__(self, other):
if isinstance(other, Atom):
return (
self.a == other.a
and self.c == other.c
and self.alpha == other.alpha
and self.beta == other.beta
and self.gamma == other.gamma
)
else:
return False
# checks if two class instances are not equal w.r.t. name, radius & valence
def __neq__(self, other):
return not self.__eq__(other)
# reprinting of class object
def __repr__(self):
return str("Unitcell(a=%s, c=%s, alpha=%s, beta=%s, gamma=%s, basis=%s)") % (
self.a,
self.c,
self.alpha,
self.beta,
self.gamma,
self.basis,
)
# redefine hash method for checking uniqueness of class instance
def __hash__(self):
return hash(self.__repr__())
# set basis vectors from (3x3) matrix in fractional coordinates
[docs]
def set_vectors(self, m3x3):
self.basis = m3x3
# set a lattice parameter
[docs]
def set_a(self, a):
"""
Description
-----------
Set the magnitude of the in-plane lattice vector a in Angstroms.
Parameters
----------
a: float
The magnitude of the in-plane lattice vector in Angstroms
Notes
-----
To retrieve a in terms of Angstroms use 'unitcell.a', whereas the
internal parameter 'unitcell._a' converts a into Bohr radii
(1 Bohr = 0.529Å), which is used for the muffin-tin potential
calculations in libphsh (CAVPOT subroutine).
"""
self.a = float(a)
self._a = self.a / 0.529 # (1 Bohr = 0.529Å)
# set c lattice parameter
[docs]
def set_c(self, c):
"""
Description
-----------
Set the magnitude of the out-of-plane lattice vector a.
Parameters
----------
c : float
The magnitude of the in-plane lattice vector in Angstroms
Notes
-----
To retrieve c in terms of Angstroms use 'unitcell.c', whereas the
internal parameter 'unitcell._c' converts c into Bohr radii
(1 Bohr = 0.529Å), which is used for the muffin-tin potential
calculations in libphsh (CAVPOT subroutine).
"""
self.c = float(c)
self._c = self.c / 0.529 # (1 Bohr = 0.529Å)
# set angle alpha in degrees
[docs]
def set_alpha(self, alpha):
self.alpha = float(alpha) % 360.0
# set angle beta in degrees
[docs]
def set_beta(self, beta):
self.beta = float(beta) % 360.0
# set angle gamma in degrees
[docs]
def set_gamma(self, gamma):
self.gamma = float(gamma) % 360.0
[docs]
class CoordinatesError(Exception):
"""Coordinate exception to raise and log duplicate coordinates."""
def __init__(self, msg):
super(CoordinatesError).__init__(type(self))
self.msg = "CoordinatesError: %s" % msg
def __str__(self):
return self.msg
def __unicode__(self):
return self.msg
[docs]
class Model(object):
"""
Generic model class.
"""
def __init__(self, unitcell, atoms, **kwargs):
"""
Constructor for Model class.
Parameters
----------
unitcell : Unitcell
An instance of the Unitcell class.
atoms : list
Array of Atom class instances which constitute the model.
"""
self.atoms = []
self.set_atoms(atoms)
self.set_unitcell(unitcell)
self.__dict__.update(kwargs)
# checks if two models are equal
def __eq__(self, other):
if isinstance(other, Atom):
return self.atoms == other.atoms and self.unitcell == other.unitcell
else:
return False
# checks if two models are not equal
def __neq__(self, other):
return not self.__eq__(other)
# reprinting of Atom object
def __repr__(self):
return str("Model(atoms=%s, unitcell=%s)") % (self.atoms, self.unitcell)
# redefine hash method for checking uniqueness of class instance
def __hash__(self):
return hash(self.__repr__())
# estimate number of inequivalent atoms
[docs]
def _nineq_atoms(self):
"""
Description
-----------
Internal method for estimating the number of inequivalent atoms
Returns
-------
nineq_atoms, element_dict : tuple
nineq_atoms : The estimated number of inequivalent atoms based on
the valence and radius of each atom.
element_dict : a dictionary of each element in the atom list where
each element contains an atom dictionary of 'nineq_atoms',
'n_atoms' and a complete 'atom_list'
Example
-------
>>> C1 = Atom('C', [0, 0, 0])
>>> Re1 = Atom('Re', [0, 0, 0], valence=2.0)
>>> Re2 = Atom('Re', [0, 0, 0], radius=1)
>>> uc = Unitcell(1, 2, [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
>>> mtz = MTZ_model(uc, atoms=[C1, Re1, Re2])
>>> print(mtz._nineq_atoms())
(3, {'Carbon': {'n_atoms': 1, 'atom_list': [Atom(Carbon, tag='C',
coordinates=[0, 0, 0], Z=6, radius=0.91, valence=0)], 'nineq_atoms':
1,'nineq_atoms_list': set([Atom(Carbon, tag='C',
coordinates=[0, 0, 0], Z=6, radius=0.91, valence=0)])}, 'Rhenium': {
'n_atoms': 2, 'atom_list': [Atom(Rhenium, tag='Re',
coordinates=[0, 0, 0], Z=75, radius=1.97, valence=2.0), Atom(Rhenium,
tag='Re', coordinates=[0, 0, 0], Z=75, radius=1, valence=0)],
'nineq_atoms': 2, 'nineq_atoms_list': set([Atom(Rhenium, tag='Re',
coordinates=[0, 0, 0], Z=75, radius=1.97, valence=2.0),
Atom(Rhenium, tag='Re', coordinates=[0, 0, 0], Z=75, radius=1,
valence=0)])}})
"""
nineq_atoms = 0
element_dict = {}
atom_dict = {}
# loop through atom list, testing each element for duplicates
# get list of elements
elements = set([atom.name for atom in self.atoms])
for element in elements:
atoms = [atom for atom in self.atoms if atom.name == element]
n_atoms = len(set(atoms))
nineq_atoms += n_atoms
atom_dict = {
"nineq_atoms": n_atoms,
"n_atoms": len(atoms),
"atom_list": atoms,
}
element_dict[element] = atom_dict
return nineq_atoms, element_dict
[docs]
def add_atom(self, element, position, **kwargs):
"""
Append an Atom instance to the model
Parameters
----------
element : str or int
Either an element name, symbol or atomic number.
position : list(float) or ndarray
(1x3) array of the fractional coordinates of the atom
within the unit cell in terms of the lattice vector a.
"""
self.atoms.append(Atom(element, position, kwargs))
[docs]
def check_coordinates(self):
"""
Check for duplicate coordinates of different atoms in model.
Raises
------
CoordinateError : exception
If duplicate positions found.
"""
positions = [str(atom.coordinates) for atom in self.atoms]
info = ""
for position in {position for position in positions if positions.count(position) > 1}:
for i, atom in enumerate([atom for atom in self.atoms if str(atom.coordinates) == position]):
info += "%s, coordinates=%s, index=%i\n" % (
str(atom),
atom.coordinates,
i,
)
if len(set(positions)) < len(self.atoms):
raise CoordinatesError("Not every atom position in model is unique!\n%s\n" % info)
@property
def name(self):
"""Returns the model name"""
return getattr(self, "_name", "")
@name.setter
def name(self, name):
"""Sets the name of the model"""
self._name = str(name) if not isinstance(name, str) else name
@property
def elements(self):
"""Returns a list of unique elements within the model"""
return {atom.element for atom in self.atoms}
@property
def atoms(self):
"""Returns a list of atoms within this model"""
if not hasattr(self, "_atoms"):
self._atoms = []
return self._atoms
@atoms.setter
def atoms(self, atoms):
"""
Set the atoms for the model.
Parameters
----------
atoms : list(Atoms)
Array of Atom instances. Entries in the list which are
not of type Atom will be ignored.
Raises
------
TypeError : exception
If atoms parameter is not a list.
"""
if isinstance(atoms, list):
self._atoms = [atom for atom in atoms if isinstance(atom, Atom)]
else:
raise TypeError
[docs]
def set_atoms(self, atoms):
"""Set atoms via the property setter for backwards compatibility."""
self.atoms = atoms
[docs]
def set_unitcell(self, unitcell):
"""
Set the unitcell for the model
Parameters
----------
unitcell : Unitcell
Instance of Unitcell class to set to model.
Raises
------
TypeError : exception
If unitcell parameter is not a Unitcell.
"""
if isinstance(unitcell, Unitcell):
self.unitcell = unitcell
else:
raise TypeError
[docs]
class MTZ_model(Model):
"""
Muffin-tin potential Model subclass for producing input file
for muffin-tin calculations in the Barbieri/Van Hove phase
shift calculation package.
"""
def __init__(self, unitcell, atoms, **kwargs):
"""
Constructor for Model class.
Parameters
----------
unitcell : Unitcell
An instance of the Unitcell class.
atoms : list
Array of Atom class instances which constitute the model.
nh : int, optional
Parameter for estimating the muffin-tin zero (default: 10).
exchange : float, optional
Hartree type exchange term alpha (default: 0.7200).
c : float, optional
The height of the slab in Angstroms - if this is much larger
than the bulk c distance then there will be a large vacuum
and therefore should be used when calculating a thin slab
rather than a bulk muffin-tin potential. Default is to lookup
the out-of-plane basis vector bulk value.
nform : int or str, optional
The phase shift calculation type, which can be 0 or 'cav' for
using the cavpot subroutine, 1 or 'wil' for using the William's
method, and 2 or 'rel' for using relativistic calculations suitable
for the CLEED package.
"""
self.atoms = []
self.set_atoms(atoms)
self.set_unitcell(unitcell)
self.set_exchange(0.72)
self.set_nh(10)
self.mtz = None
self.__dict__.update(kwargs)
[docs]
def set_nh(self, nh):
"""Sets the nh muffin-tin zero estimation parameter"""
self.nh = int(nh) # check this is not float
[docs]
def set_exchange(self, alpha):
"""Sets the alpha exchange term for muffin-tin calculation"""
try:
self.exchange = float(alpha)
except (TypeError, ValueError):
return
# set form of muffin-tin calculation: 0=cav, 1=wil, 2=rel
[docs]
def set_slab_c(self, c):
"""Set the maximum height of the slab in Angstroms.
If this is much larger than the bulk c distance then there will be a
large vacuum and therefore should be used when calculating a thin slab
rather than a bulk muffin-tin potential.
Examples
--------
For Re the bulk c distance is 2.76Å, whereas a possible slab c distance
could be ~10Å.
"""
try:
self.c = float(c)
self._c = self.c / 0.529
except (TypeError, ValueError):
return
[docs]
def load_from_file(self, filename):
"""
Description
-----------
Load an input file and update the class instance variables
Parameters
----------
filename : str
The path of the input file (e.g. cluster*.i or *slab*.i)
Raises
------
IOError : exception
If the file cannot be read.
TypeError : exception
If a input line cannot be parsed correctly.
"""
from phaseshifts.leed import Converter, CLEEDInputValidator
filename = glob(os.path.expanduser(os.path.expandvars(filename)))[0]
if CLEEDInputValidator.is_cleed_file(filename):
cleed_model = Converter.import_CLEED(filename)
self.__dict__.update(cleed_model.__dict__)
return
self._load_input_file(filename)
[docs]
def create_atorbs(self, **kwargs):
"""
Description
-----------
Create Atorb input files for each element in model.
Arguments
---------
output_dir : str
Path to output directory for files
library_dir : str
Path to library of input files. Use this if you've previously
created input files for a given element as it doesn't need to
recalculate the radial charge density every time - note this
is also a workaround to allow precalculated files if your machine's
output is drastically different from what is expected.
config : dict
See help(atorb.gen_input) for list of keywords and values.
Returns
-------
output_files : dict
Dictionary list of atorb*.i input files for the Atorb class to
calculate the charge density from.
"""
# check output path
if "output_dir" in kwargs:
output_dir = kwargs["output_dir"]
else:
if "library_dir" in kwargs:
output_dir = kwargs["library_dir"]
else:
output_dir = os.path.abspath(".")
atorb_files = []
at_files = []
# generate input files for atomic orbitals and charge densities
for element in set([str(atom.element.symbol) for atom in self.atoms]):
if not os.path.isdir(os.path.join(output_dir, "Atorb")):
try:
os.makedirs(os.path.join(output_dir, "Atorb"))
except OSError as e:
# Only ignore "directory already exists" errors, re-raise others
if e.errno != getattr(errno, "EEXIST", 17):
raise
atorb_file = os.path.join(output_dir, "Atorb", "atorb_%s.i" % element)
if not os.path.isfile(atorb_file): # create new atorb input file
atorb_file = atorb.Atorb.gen_input(element, filename=atorb_file)
at_file = os.path.join(output_dir, "Atorb", "at_" + element + ".i")
if not os.path.isfile(at_file):
at_file = atorb.Atorb.calculate_Q_density(
input=atorb_file, output_dir=os.path.join(output_dir, "Atorb")
)
# update lists
atorb_files.append(atorb_file)
at_files.append(at_file)
return {"atorb_files": atorb_files, "at_files": at_files}
[docs]
def gen_atomic(self, **kwargs):
"""
Description
-----------
Create atomic*.i input file for MTZ input based on model or
list of files.
Parameters
----------
input_dir : str
Input directory where at*.i files are kept.
input_files : tuple
List of input files to generate atomic input file from.
output_file : str
The filename of the resulting atomic*.i output file, which is
simply a superimposed set of the radial charge densities from
the individual input files.
Returns
-------
output_file : str
Returns the output file path string.
Raises
------
IOError : exception
If either input directory or files do not exist.
Notes
-----
If 'input_files' is not given then the default list of input files are
inferred from the list of atoms in the model.
"""
# input
if "input_dir" in kwargs:
input_dir = os.path.abspath(glob(os.path.expanduser(os.path.expandvars(kwargs["input_dir"])))[0])
else:
input_dir = os.path.abspath(".")
if os.path.isfile(input_dir):
input_dir = os.path.dirname(input_dir)
if not os.path.isdir(input_dir):
raise IOError("'%s' is not a valid input directory - " "does not exist!" % input_dir)
# output filename
if "output_file" in kwargs:
output_file = os.path.abspath(kwargs["output_file"])
else:
output_file = os.path.abspath("atomic.i")
# get list of input
if "input_files" in kwargs:
files = kwargs["input_files"]
else: # assume using atoms from model
files = [os.path.join(input_dir, "at_" + atom.element.symbol + ".i") for atom in self.atoms]
# generate atomic.i input file by appending multiple at.i files
with open(output_file, "w") as f:
# loop through each atomic charge density file in list
for input_file in files:
if not os.path.isfile(str(input_file)) or input_file is None:
raise IOError("Radial charge density file " "'%s' does not exist!" % input_file)
# append next input file to output
with open(input_file) as infile:
f.write(infile.read())
return output_file
[docs]
def get_MTZ(self, filename):
"""Retrieves muffin-tin potential from file"""
try:
with open(filename, "r") as f:
self.mtz = float([line for line in f][0]) # read first line
except IOError:
raise IOError
[docs]
def calculate_MTZ(self, mtz_string="", **kwargs):
"""
Description
-----------
Calculate the muffin-tin potential (MTZ) for a given cluster file
Parameters
----------
atomic_file : str
The path to the atomic input file. If this is omitted the default
is generate one using the MTZ_model.gen_atomic() method.
cluster_file : str
The path to the cluster input file which may be a bulk or slab
model.
slab : int or bool
Determines whether the MTZ calculation is for a slab model (True).
The default is a bulk calculation.
output : dict
Dictionary output of 'mtz' - muffin-tin potential & 'output_file'
- the path to the MTZ output file.
Returns
-------
output_files : list(str)
Paths to the MTZ output file.
"""
# check to see if cluster input exists
if "cluster_file" in kwargs:
cluster_file = os.path.abspath(glob(os.path.expanduser(os.path.expandvars(kwargs["cluster_file"])))[0])
else:
cluster_file = os.path.abspath("cluster.i")
if not os.path.isfile(cluster_file):
raise IOError("MTZ cluster file '%s' does not exist!" % cluster_file)
if not os.access(os.path.dirname(cluster_file), os.W_OK) and "atomic_file" not in kwargs:
raise IOError("Do not have write access to '%s'" % os.path.dirname(cluster_file))
# determine type of calculation - bulk or slab
if "slab" in kwargs:
slab = bool(kwargs["slab"])
fid = "mtz"
else:
slab = False
fid = "bmtz"
mtz_string = ""
if "atomic_file" in kwargs:
atomic_file = os.path.abspath(kwargs["atomic_file"])
if not os.path.isfile(atomic_file):
raise IOError("Appended radial charge densities file " "'%s' does not exist!" % atomic_file)
else: # generate on the fly
input_dir = os.path.abspath(os.path.dirname(cluster_file))
self.create_atorbs(output_dir=input_dir)
self.gen_atomic(input_dir=input_dir)
if "output_file" in kwargs:
output_file = os.path.abspath(kwargs["output_file"])
else:
output_file = "%s.i" % fid
try:
os.makedirs(os.path.dirname(output_file))
except WindowsError:
pass
if os.path.isfile(output_file):
move(output_file, output_file + ".bak")
# create mufftin debug file
mufftin_file = os.path.splitext(output_file)[0] + "_mufftin.d"
info_file = os.path.splitext(output_file)[0] + "_info.txt"
# call cavpot routine
self.mtz = libphsh.cavpot(
mtz_string,
int(slab),
atomic_file,
cluster_file,
mufftin_file,
output_file,
info_file,
)
# check to see if new file has been written
if not os.path.isfile(output_file):
raise IOError("Failed to write muffin-tin potential file '%s'" % output_file)
return output_file
[docs]
def get_elements(self):
"""Return the unique elements in model"""
return set([atom.name for atom in self.atoms])
# #==============================================================================
# # Testing
# #==============================================================================
# at = Atom('C', [0, 0, 0])
# ab = Atom('Re', [0, 0, 0], tag='Re2')
# ac = Atom('Re', [0, 0, 0], tag='Re1')
# uc = Unitcell(1, 2, [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
# print(uc)
# print(set([at, ab, ac]))
#
# mtz = MTZ_model(uc, atoms=[at, ab, ac])
# mtz.load_from_file(
# 'C:\\Users\\Liam\\Dropbox\\Programming\\Python\\LEED-PyV\\phaseshifts\\test\\Re0001\\cluster_Re_bulk.i'
# )
# print(mtz.get_elements())
# mtz.load_from_file('C:\\Users\\kss07698\\Desktop\\test_cluster.bak.i')
# mtz.gen_input(filename='C:\\Users\\kss07698\\Desktop\\test_cluster.bak.i')