#!/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 elements
import atorb
import os
from copy import deepcopy
from shutil import move
from glob import glob
from 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.], **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):
try:
self.coordinates = coordinates
self._coordinates = [r / 0.529 for r in coordinates]
except any as e:
raise e
# 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:
pass
[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 set([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)
[docs] def set_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_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:
pass
# set form of muffin-tin calculation: 0=cav, 1=wil, 2=rel
[docs] def set_slab_c(self, c):
"""
Description
-----------
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:
pass
[docs] def load_from_file(self, filename):
"""
Description
-----------
Load an cluster/slab 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.
"""
filename = glob(os.path.expanduser(os.path.expandvars(filename)))[0]
try:
with open(filename, 'r') as f:
self.header = f.readline()
a = float(f.readline().split('#')[0].split()[0]) * 0.529
a1 = [t(s) for (t, s) in zip((float, float, float),
f.readline().split('#')[0].split()[:3])]
a2 = [t(s) for (t, s) in zip((float, float, float),
f.readline().split('#')[0].split()[:3])]
a3 = [t(s) for (t, s) in zip((float, float, float),
f.readline().split('#')[0].split()[:3])]
basis = [a1, a2, a3]
c = float(a3[-1]) * 0.529 # change to Angstroms from Bohr
self.set_unitcell(Unitcell(a, c, basis))
self.nineq_atoms = int(f.readline().split('#')[0].split()[0])
self.atoms = []
line = f.readline().split('#')[0]
while line.split()[0].isalpha():
n, Z, val, rad = [t(s) for (t, s) in zip(
(int, float, float, float),
f.readline().split('#')[0].split()[:4])]
for i in range(0, n):
position = [t(s) for (t, s) in
zip((float, float, float),
f.readline().split('#')[0].split()[:3])]
atom = Atom(line.split()[0].capitalize(),
coordinates=position,
tag=line.split()[1], valence=val,
radius=rad, Z=int(Z))
self.atoms.append(atom)
line = f.readline().split('#')[0]
#print(line)
self.set_nform(line.split()[0])
self.set_exchange(f.readline().split('#')[0].split()[0])
self.set_nh(f.readline().split('#')[0].split()[0])
except IOError:
raise IOError("cannot read '%s'" % filename)
except ValueError:
raise ValueError("malformatted input in '%s'" % 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 WindowsError:
pass # directory already exists on Windows
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 == 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')