##################################################################
## (c) Copyright 2015- by Jaron T. Krogel ##
##################################################################
#====================================================================#
# structure.py #
# Support for atomic structure I/O, generation, and manipulation. #
# #
# Content summary: #
# Structure #
# Represents a simulation cell containing a set of atoms. #
# Many functions for manipulating structures or obtaining #
# data regarding local atomic structure. #
# #
# generate_cell #
# User-facing function to generate an empty simulation cell. #
# #
# generate_structure #
# User-facing function to specify arbitrary atomic structures #
# or generate structures corresponding to atoms, dimers, or #
# crystals. #
# #
#====================================================================#
"""
The :py:mod:`structure` module provides support for atomic structure I/O,
generation, and manipulation.
List of module contents
-----------------------
Read cif file functions:
* :py:func:`read_cif_celldata`
* :py:func:`read_cif_cell`
* :py:func:`read_cif`
Operations on logical conditions:
* :py:func:`equate`
* :py:func:`negate`
Create a Monkhorst-Pack k-point mesh function
* :py:func:`kmesh`
Tile matrix malipulation functions
* :py:func:`reduce_tilematrix`
* :py:func:`tile_magnetization`
Rotate plane function
* :py:func:`rotate_plane`
Trivial filter function
* :py:func:`trivial_filter`
* :py:class:`MaskFilter`
* :py:func:`optimal_tilematrix`
Base class for :py:class:`Structure` class:
* :py:class:`Sobj`
Base class for :py:class:`DefectStructure`, :py:class:`Crystal`, and :py:class:`Jellium` classes:
* :py:class:`Structure`
SeeK-path functions
* :py:func:`_getseekpath`
* :py:func:`get_conventional_cell`
* :py:func:`get_primitive_cell`
* :py:func:`get_kpath`
* :py:func:`get_symmetry`
* :py:func:`get_structure_with_bands`
* :py:func:`get_band_tiling`
* :py:func:`get_seekpath_full`
Interpolate structures functions
* :py:func:`interpolate_structures`
Animate structures functions
* :py:func:`structure_animation`
Concrete :py:class:`Structure` classes:
* :py:class:`DefectStructure`
* :py:class:`Crystal`
* :py:class:`Jellium`
Structure generation functions:
* :py:func:`generate_cell`
* :py:func:`generate_structure`
* :py:func:`generate_atom_structure`
* :py:func:`generate_dimer_structure`
* :py:func:`generate_trimer_structure`
* :py:func:`generate_jellium_structure`
* :py:func:`generate_crystal_structure`
* :py:func:`generate_defect_structure`
Read structure functions
* :py:func:`read_structure`
Module contents
---------------
"""
import os
from pathlib import Path
import numpy as np
from copy import deepcopy
from random import randint
import itertools
from numpy import (
cos,
cross,
dot,
pi,
sin,
sqrt,
)
from numpy.linalg import inv, det, norm
import numpy.typing as npt
from .unit_converter import convert
from .numerics import nearest_neighbors, convex_hull, voronoi_neighbors
from .periodic_table import Elements
from .fileio import XsfFile, PoscarFile
from .developer import DevBase, obj, unavailable, error
from .utilities import path_string
from . import numpy_extensions as npe
try:
from scipy.special import erfc
except:
erfc = unavailable('scipy.special','erfc')
#end try
try:
import matplotlib.pyplot as plt
from matplotlib.pyplot import plot,subplot,title,xlabel,ylabel
except:
plot,subplot,title,xlabel,ylabel,plt = unavailable('matplotlib.pyplot','plot','subplot','title','xlabel','ylabel','plt')
#end try
# installation instructions for spglib interface
#
# this is bootstrapped off of spglib's ASE Python interface
#
# installation of spglib
# go to http://sourceforge.net/projects/spglib/files/
# click on Download spglib-1.8.2.tar.gz (952.6 kB)
# unpack directory (tar -xzf spglib-1.8.2.tar.gz)
# enter ase directory (cd spglib-1.8.2/python/ase/)
# build and install (sudo python setup.py install)
#try:
# from pyspglib import spglib
#except:
# spglib = unavailable('pyspglib','spglib')
##end try
try:
import spglib
except:
spglib = unavailable('spglib')
#end try
# installation instructions to enable cif file read
#
# cif file support in Nexus currently requires two external libraries
# PyCifRW - base interface to read cif files into object format: CifFile
# cif2cell - translation layer from CifFile object to cell reconstruction: CellData
# (note: cif2cell installation includes PyCifRW)
#
# installation of cif2cell
# go to http://sourceforge.net/projects/cif2cell/
# click on Download (example: cif2cell-1.2.10.tar.gz)
# unpack directory (tar -xzf cif2cell-1.2.10.tar.gz)
# enter directory (cd cif2cell-1.2.10)
# install cif2cell (python setup.py install)
# check python installation
# >python
# >>>from CifFile import CifFile
# >>>from uctools import CellData
#
# Nexus is currently compatible with
# cif2cell-1.2.10 and PyCifRW-3.3
# cif2cell-1.2.7 and PyCifRW-4.1.1
# compatibility last tested: 20 Mar 2017
#
try:
from CifFile import CifFile
except:
CifFile = unavailable('CifFile','CifFile')
#end try
try:
from cif2cell.uctools import CellData
except:
CellData = unavailable('cif2cell.uctools','CellData')
#end try
cif2cell_unit_dict = dict(angstrom='A',bohr='B',nm='nm')
[docs]
def read_cif_celldata(filepath,block=None,grammar='1.1'):
# read cif file with PyCifRW
path,cif_file = os.path.split(filepath)
if path!='':
cwd = os.getcwd()
os.chdir(path)
#end if
cf = CifFile(cif_file,grammar=grammar)
#cf = ReadCif(cif_file,grammar=grammar)
if path!='':
os.chdir(cwd)
#end if
if block is None:
block = list(cf.keys())[0]
#end if
cb = cf.get(block)
if cb is None:
error('block {0} was not found in cif file {1}'.format(block,filepath),'read_cif_celldata')
#end if
# repack H-M symbols as normal strings so CellData.getFromCIF won't choke on unicode
#for k in ['_symmetry_space_group_name_H-M','_space_group_name_H-M_alt','_symmetry_space_group_name_h-m','_space_group_name_h-m_alt']:
# if k in cb.block:
# v = cb.block[k]
# if isinstance(v,(list,tuple)):
# for i in range(len(v)):
# if isinstance(v[i],unicode):
# v[i] = str(v[i])
# #end if
# #end for
# #end if
# #end if
##end for
# extract structure from CifFile with uctools CellData class
cd = CellData()
cd.getFromCIF(cb)
return cd
#end def read_cif_celldata
[docs]
def read_cif_cell(filepath,block=None,grammar='1.1',cell='prim'):
cd = read_cif_celldata(filepath,block,grammar)
if cell.startswith('prim'):
cell = cd.primitive()
elif cell.startswith('conv'):
cell = cd.conventional()
else:
error('cell argument must be primitive or conventional\nyou provided: {0}'.format(cell),'read_cif_cell')
#end if
return cell
#end def read_cif_cell
[docs]
def read_cif(filepath,block=None,grammar='1.1',cell='prim',args_only=False):
if isinstance(filepath,str):
cell = read_cif_cell(filepath,block,grammar,cell)
else:
cell = filepath
#end if
# create Structure object from cell
if cell.alloy:
error('cannot handle alloys','read_cif')
#end if
units = cif2cell_unit_dict[cell.unit]
scale = float(cell.lengthscale)
scale = convert(scale,units,'A')
units = 'A'
axes = scale*np.array(cell.latticevectors,dtype=float)
elem = []
pos = []
for wyckoff_atoms in cell.atomdata:
for atom in wyckoff_atoms:
elem.append(str(list(atom.species.keys())[0]))
pos.append(atom.position)
#end for
#end for
pos = dot(np.array(pos,dtype=float),axes)
if not args_only:
s = Structure(
axes = axes,
elem = elem,
pos = pos,
units = units
)
return s
else:
return axes,elem,pos,units
#end if
#end def read_cif
[docs]
def equate(expr):
return expr
#end def equate
[docs]
def negate(expr):
return not expr
#end def negate
[docs]
def kmesh(kaxes,dim,shift=None):
'''
Create a Monkhorst-Pack k-point mesh
'''
if shift is None:
shift = (0.,0,0)
#end if
ndim = len(dim)
d = np.array(dim)
s = np.array(shift)
npe.reshape_inplace(s, (1, ndim))
npe.reshape_inplace(d, (1, ndim))
kp = np.empty((1,ndim),dtype=float)
kgrid = np.empty((d.prod(),ndim))
n=0
for k in range(dim[2]):
for j in range(dim[1]):
for i in range(dim[0]):
kp[:] = i,j,k
kp = dot((kp+s)/d,kaxes)
#kp = (kp+s)/d
kgrid[n] = kp
n+=1
#end for
#end for
#end for
return kgrid
#end def kmesh
[docs]
def reduce_tilematrix(tiling):
tiling = np.array(tiling)
t = np.array(tiling,dtype=int)
if np.abs(tiling-t).sum()>1e-6:
Structure.class_error('requested tiling is non-integer\n tiling requested: '+str(tiling))
#end if
dim = len(t)
matrix_tiling = t.shape == (dim,dim)
if matrix_tiling:
if np.abs(det(t))==0:
Structure.class_error('requested tiling matrix is singular\ntiling requested: {0}'.format(t))
#end if
#find a tiling tuple from the tiling matrix
# do this by shearing the tiling matrix (or equivalently the tiled cell)
# until it is orthogonal (in the untiled cell axes)
# this is just rearranging triangular tiles of space to reshape the cell
# so that t1*t2*t3 = det(T) = det(A_tiled)/det(A_untiled)
#this way the atoms in the (perhaps oddly shaped) supercell can be
# obtained from simple translations of the untiled cell positions
T = t #tiling matrix
tilematrix = T.copy()
del t
tbar = np.identity(dim) #basis for shearing
dr = list(range(dim))
#dr = [1,0,2]
other = dim*[0] # other[d] = dimensions other than d
for d in dr:
other[d] = set(dr)-set([d])
#end for
#move each axis to be parallel to barred directions
# these are volume preserving shears of the supercell
# each shear keeps two cell face planes fixed while moving the others
tvecs = []
for dp in [(0,1,2),(2,0,1),(1,2,0),(2,1,0),(0,2,1),(1,0,2)]:
success = True
Tnew = np.array(T,dtype=float) #sheared/orthogonal tiling matrix
for d in dr:
tb = tbar[dp[d]]
t = T[d]
d2,d3 = other[d]
n = cross(Tnew[d2],Tnew[d3]) #vector normal to 2 cell faces
vol = dot(n,t)
bcomp = dot(n,tb)
if np.abs(bcomp)<1e-6:
success = False
break
#end if
tn = vol*1./bcomp*tb #new axis vector
Tnew[d] = tn
#end for
if success:
# apply inverse permutation, if needed
Tn = Tnew.copy()
for d in dr:
d2 = dp[d]
Tnew[d2] = Tn[d]
#end for
#the resulting tiling matrix should be diagonal and integer
tr = np.diag(Tnew)
nondiagonal = np.abs(Tnew-np.diag(tr)).sum()>1e-6
if nondiagonal:
Structure.class_error('could not find a diagonal tiling matrix for generating tiled coordinates')
#end if
tvecs.append(np.abs(tr))
#end if
#end for
tvecs_old = tvecs
tvecs = []
tvset = set()
for tv in tvecs_old:
tvk = tuple(np.array(np.around(1e7*tv),dtype=np.uint64))
if tvk not in tvset:
tvset.add(tvk)
tvecs.append(tv)
#end if
#end for
tilevector = np.array(tvecs)
else:
tilevector = t
tilematrix = np.diag(t)
#end if
return tilematrix,tilevector
#end def reduce_tilematrix
[docs]
def rotate_plane(plane,angle,points,units='degrees'):
if units=='degrees':
angle *= pi/180
elif not units.startswith('rad'):
error('angular units must be degrees or radians\nyou provided: {0}'.format(angle),'rotate_plane')
#end if
c = cos(angle)
s = sin(angle)
if plane=='xy':
R = [[ c,-s, 0],
[ s, c, 0],
[ 0, 0, 1]]
elif plane=='yx':
R = [[ c, s, 0],
[-s, c, 0],
[ 0, 0, 1]]
elif plane=='yz':
R = [[ 1, 0, 0],
[ 0, c,-s],
[ 0, s, c]]
elif plane=='zy':
R = [[ 1, 0, 0],
[ 0, c, s],
[ 0,-s, c]]
elif plane=='zx':
R = [[ c, 0, s],
[ 0, 1, 0],
[-s, 0, c]]
elif plane=='xz':
R = [[ c, 0,-s],
[ 0, 1, 0],
[ s, 0, c]]
else:
error('plane must be xy/yx/yz/zy/zx/xz\nyou provided: {0}'.format(plane),'rotate_plane')
#end if
R = np.array(R,dtype=float)
return dot(R,points.T).T
#end def rotate_plane
[docs]
def recenter_points(pos, center, axes):
"""Center a given set of points in a cell around a provided center.
Parameters
----------
pos : NDArray
Array of *N* positions with shape (*N*,3)
center : NDArray
Position of the center
axes : NDArray
Array of the cell vectors with shape (3,3)
Returns
-------
pos : NDArray
Array of positions centered around the given center
Notes
-----
This function also ensures that points close (within 1e-12) to the minimum edge (-0.5) of the
cell are placed exactly on that edge. The intent here is to make sure that atoms close to or
on the leading edge (+0.5) are wrapped around to retain periodicity.
"""
axinv = inv(axes)
for i in range(len(pos)):
u = dot(pos[i] - center, axinv)
u -= np.floor(u+.5)
u[np.abs(u-.5)<1e-12] = -0.5
u[np.abs(u+.5)<1e-12] = -0.5
pos[i] = dot(u,axes) + center
#end for
return pos
opt_tm_matrices = obj()
opt_tm_wig_indices = obj()
[docs]
def trivial_filter(T):
return True
#end def trival_filter
[docs]
class MaskFilter(DevBase):
[docs]
def set(self,mask,dim=3):
omask = np.array(mask)
mask = np.array(mask,dtype=bool)
if mask.size==dim:
mvec = mask.ravel()
mask = np.empty((dim,dim),dtype=bool)
i=0
for mi in mvec:
j=0
for mj in mvec:
mask[i,j] = mi==mj
j+=1
#end for
i+=1
#end for
elif mask.shape!=(dim,dim):
error(
"shape of mask array must be {0},{0}\n"
"shape received: {1},{2}\n"
"mask array received: {3}".format(
dim,mask.shape[0],mask.shape[1],omask),'optimal_tilematrix'
)
#end if
self.mask = mask==False
#end def set
def __call__(self,T):
return (T[self.mask]==0).all()
#end def __call__
#end class MaskFilter
mask_filter = MaskFilter()
[docs]
def optimal_tilematrix(axes,volfac,dn=1,tol=1e-3,filter=trivial_filter,mask=None,nc=5,Tref=None):
if mask is not None:
mask_filter.set(mask)
filter = mask_filter
#end if
dim = 3
if isinstance(axes,Structure):
axes = axes.axes
else:
axes = np.array(axes,dtype=float)
#end if
if not isinstance(volfac,int):
volfac = int(np.around(volfac))
#end if
volume = np.abs(det(axes))*volfac
axinv = inv(axes)
cube = volume**(1./3)*np.identity(dim)
if Tref is None:
Tref = np.array(np.around(dot(cube,axinv)),dtype=int)
else:
Tref = np.asarray(Tref)
#end if
# calculate and store all tiling matrix variations
if dn not in opt_tm_matrices:
mats = []
rng = tuple(range(-dn,dn+1))
for n1 in rng:
for n2 in rng:
for n3 in rng:
for n4 in rng:
for n5 in rng:
for n6 in rng:
for n7 in rng:
for n8 in rng:
for n9 in rng:
mats.append((n1,n2,n3,n4,n5,n6,n7,n8,n9))
#end for
#end for
#end for
#end for
#end for
#end for
#end for
#end for
#end for
mats = np.array(mats,dtype=int)
npe.reshape_inplace(mats, ((2*dn+1)**(dim*dim), dim, dim))
opt_tm_matrices[dn] = mats
else:
mats = opt_tm_matrices[dn]
#end if
# calculate and store all wigner image indices
if nc not in opt_tm_wig_indices:
inds = []
rng = tuple(range(-nc,nc+1))
for k in rng:
for j in rng:
for i in rng:
if i!=0 or j!=0 or k!=0:
inds.append((i,j,k))
#end if
#end for
#end for
#end for
inds = np.array(inds,dtype=int)
opt_tm_wig_indices[nc] = inds
else:
inds = opt_tm_wig_indices[nc]
#end if
# track counts of tiling matrices
ntilings = len(mats)
nequiv_volume = 0
nfilter = 0
nequiv_inscribe = 0
nequiv_wigner = 0
nequiv_cubicity = 0
nequiv_shape = 0
# try a faster search for cells w/ target volume
det_inds_p = [
[(0,0),(1,1),(2,2)],
[(0,1),(1,2),(2,0)],
[(0,2),(1,0),(2,1)]
]
det_inds_m = [
[(0,0),(1,2),(2,1)],
[(0,1),(1,0),(2,2)],
[(0,2),(1,1),(2,0)]
]
volfacs = np.zeros((len(mats),),dtype=int)
for (i1,j1),(i2,j2),(i3,j3) in det_inds_p:
volfacs += (Tref[i1,j1]+mats[:,i1,j1])*(Tref[i2,j2]+mats[:,i2,j2])*(Tref[i3,j3]+mats[:,i3,j3])
#end for
for (i1,j1),(i2,j2),(i3,j3) in det_inds_m:
volfacs -= (Tref[i1,j1]+mats[:,i1,j1])*(Tref[i2,j2]+mats[:,i2,j2])*(Tref[i3,j3]+mats[:,i3,j3])
#end for
Tmats = mats[np.abs(volfacs)==volfac]
nequiv_volume = len(Tmats)
# find the set of cells with maximal inscribing radius
inscribe_tilings = []
rmax = -1e99
for mat in Tmats:
T = Tref + mat
if filter(T):
nfilter+=1
Taxes = dot(T,axes)
rc1 = norm(cross(Taxes[0],Taxes[1]))
rc2 = norm(cross(Taxes[1],Taxes[2]))
rc3 = norm(cross(Taxes[2],Taxes[0]))
r = 0.5*volume/max(rc1,rc2,rc3) # inscribing radius
if r>rmax or np.abs(r-rmax)<tol:
inscribe_tilings.append((r,T,Taxes))
rmax = r
#end if
#end if
#end for
# find the set of cells w/ maximal wigner radius out of the inscribing set
wigner_tilings = []
rwmax = -1e99
for r,T,Taxes in inscribe_tilings:
if np.abs(r-rmax)<tol:
nequiv_inscribe+=1
rw = 1e99
for ind in inds:
rw = min(rw,0.5*norm(dot(ind,Taxes)))
#end for
if rw>rwmax or np.abs(rw-rwmax)<tol:
wigner_tilings.append((rw,T,Taxes))
rwmax = rw
#end if
#end if
#end for
# find the set of cells w/ maximal cubicity
# (minimum cube_deviation)
cube_tilings = []
cmin = 1e99
for rw,T,Ta in wigner_tilings:
if np.abs(rw-rwmax)<tol:
nequiv_wigner+=1
dc = volume**(1./3)*sqrt(2.)
d1 = np.abs(norm(Ta[0]+Ta[1])-dc)
d2 = np.abs(norm(Ta[1]+Ta[2])-dc)
d3 = np.abs(norm(Ta[2]+Ta[0])-dc)
d4 = np.abs(norm(Ta[0]-Ta[1])-dc)
d5 = np.abs(norm(Ta[1]-Ta[2])-dc)
d6 = np.abs(norm(Ta[2]-Ta[0])-dc)
cube_dev = (d1+d2+d3+d4+d5+d6)/(6*dc)
if cube_dev<cmin or np.abs(cube_dev-cmin)<tol:
cube_tilings.append((cube_dev,rw,T,Ta))
cmin = cube_dev
#end if
#end if
#end for
# prioritize selection by "shapeliness" of tiling matrix
# prioritize positive diagonal elements
# penalize off-diagonal elements
# penalize negative off-diagonal elements
shapely_tilings = []
smax = -1e99
for cd,rw,T,Taxes in cube_tilings:
if np.abs(cd-cmin)<tol:
nequiv_cubicity+=1
d = np.diag(T)
o = (T-np.diag(d)).ravel()
s = np.sign(d).sum()-(np.abs(o)>0).sum()-(o<0).sum()
if s>smax or np.abs(s-smax)<tol:
shapely_tilings.append((s,rw,T,Taxes))
smax = s
#end if
#end if
#end for
# prioritize selection by symmetry of tiling matrix
ropt = -1e99
Topt = None
Taxopt = None
diagonal = []
symmetric = []
antisymmetric = []
other = []
for s,rw,T,Taxes in shapely_tilings:
if np.abs(s-smax)<tol:
nequiv_shape+=1
Td = np.diag(np.diag(T))
if np.abs(Td-T).sum()==0:
diagonal.append((rw,T,Taxes))
elif np.abs(T.T-T).sum()==0:
symmetric.append((rw,T,Taxes))
elif np.abs(T.T+T-2*Td).sum()==0:
antisymmetric.append((rw,T,Taxes))
else:
other.append((rw,T,Taxes))
#end if
#end if
#end for
s = 1
if len(diagonal)>0:
cells = diagonal
elif len(symmetric)>0:
cells = symmetric
elif len(antisymmetric)>0:
cells = antisymmetric
s = -1
elif len(other)>0:
cells = other
else:
cells = []
#end if
skew_min = 1e99
if len(cells)>0:
for rw,T,Taxes in cells:
Td = np.diag(np.diag(T))
skew = np.abs(T.T-s*T-(1-s)*Td).sum()
if skew<skew_min:
ropt = rw
Topt = T
Taxopt = Taxes
skew_min = skew
#end if
#end for
#end if
if Taxopt is None:
error(
"optimal tilematrix for volfac={0} not found with tolerance {1}\n"
"difference range (dn): {2}\n"
"tiling matrices searched: {3}\n"
"cells with target volume: {4}\n"
"cells that passed the filter: {5}\n"
"cells with equivalent inscribing radius: {6}\n"
"cells with equivalent wigner radius: {7}\n"
"cells with equivalent cubicity: {8}\n"
"matrices with equivalent shapeliness: {9}\n"
"please try again with dn={10}".format(
volfac, tol, dn, ntilings, nequiv_volume, nfilter, nequiv_inscribe,
nequiv_wigner, nequiv_cubicity, nequiv_shape, dn+1
)
)
#end if
if det(Taxopt)<0:
Topt = -Topt
#end if
return Topt,ropt
#end def optimal_tilematrix
[docs]
class Sobj(DevBase):
None
#end class Sobj
[docs]
class Structure(Sobj):
"""General class for all physical structures
Attributes
----------
axes : NDArray
Lattice vectors of the cell as a square matrix of dimension
``dim x dim``.
scale : int or float
Scaling for all other physical values.
elem : NDArray of str
Array of atomic symbols.
pos : NDArray of float
Positions of the atoms.
center : NDArray of float
Defined center of the cell.
kpoints : NDArray of float
Array containing the positions of k-points.
kweights : NDArray of float
Weight of individual k-points.
units : str
Units of the atom positions.
dim : int, default=3
Dimensionality of the Structure.
operations : Mapping of str to function
Operations to perform on the structure.
See :py:meth:`~.Structure.set_operations()` for more information.
background_charge : int, default=0
The total background charge of the system. Positive for cations,
negative for anions, and zero for neutral systems.
bconds : NDArray of str
Boundary conditions either in all directions or specified for
each dimension
mag : NDArray of float or None
Magnetic moments of each atom, or ``None`` if system is not
magnetized.
tmatrix : NDArray of int or None
A vector of ints for tiling the cell in each dimension.
See :py:meth:`~.Structure.tile()` for more information.
frozen : NDArray of bool or None, default=None
Mask array of booleans with the same length as the number of
atoms and width of ``dim``, where ``True`` indicates the atoms
are locked in place and ``False`` indicates they are free to
move.
Parameters
----------
axes : ArrayLike, optional
Lattice vectors of the cell
scale : int or float, default=1.
Scaling for all other physical values.
See :py:meth:`~.Structure.rescale()` for more information.
elem : ArrayLike, optional
Array of atomic symbols.
pos : ArrayLike, optional
Positions of the atoms.
elem_pos : str, optional
Multiline string with each line containing an atom where an atom
is an element and its position. Overrides ``elem`` and ``pos``.
mag : ArrayLike, optional
Magnetic moments of each atom.
center : ArrayLike, optional
Defined center of the cell, defaults to the ``[0.5, 0.5, 0.5]``
point of the cell.
kpoints : ArrayLike, optional
Array containing the positions of k-points.
kweights : ArrayLike, optional
Weight of individual k-points, must be as long as the number of
k-points.
kgrid : int, optional
Number of subdivisions to create a Monkhorst-Pack k-point mesh.
See :py:func:`~.kmesh` for a more in-depth explanation.
Overrides ``kpoints`` if specified.
kshift : ArrayLike, optional
Vector to translate k-points if k-grid is specified.
permute : ArrayLike, optional
Vector to permute the structure.
See :py:meth:`~.Structure.permute()` for more information.
units : str, optional
Units of the positions. See :py:mod:`unit_converter.py` for a
full list of supported units.
tiling : ArrayLike of int, optional
A vector of ints for tiling the cell in each dimension.
See :py:meth:`~.Structure.tile()` for more information.
rescale : bool, default=True
``True`` will rescale the supplied structural information with
the scaling factor provided (see ``scale``), ``False`` sets the
``scale`` without altering the provided structural information.
dim : int, default=3
Dimensionality of the Structure. See Notes for more information.
operations : iterable of str, optional
Operations to perform on the structure.
See :py:meth:`~.Structure.set_operations()` for more information.
background_charge : int, default=0
The total background charge of the system. Positive for cations,
negative for anions, and zero for neutral systems.
frozen : ArrayLike of bool, optional
Mask array of booleans with the same length as the number of
atoms and width of ``dim``.
See :py:meth:`~.Structure.set_frozen()` for more information.
bconds : str or tuple of str, optional
Boundary conditions either in all directions or specified for
each dimension. Defaults to periodic in all directions.
posu : ArrayLike, optional
The positions of the atoms in units of the lattice parameter.
Overrides ``pos``.
use_prim : bool, optional
Option to convert the unit cell to a primitive cell.
Requires that the ``seekpath`` package is installed [0]_.
add_kpath : bool, default=False
Optionally add k-points to the primitive cell. Only used if
``use_prim`` is ``True``.
symm_kgrid : bool, default=False
Option for generating a Monkhorst-Pack k-point grid with only
symmetric k-points. Requires the ``spglib`` package [1]_.
See :py:meth:`~.Structure.add_symmetrized_mesh()` for more
information.
Notes
-----
Currently the :py:class:`~.Structure` class only partially supports
2-dimensional structures.
References
----------
.. [0] https://seekpath.readthedocs.io/en/latest/
.. [1] https://spglib.readthedocs.io/en/stable/
"""
operations = obj()
[docs]
@classmethod
def set_operations(cls):
cls.operations.set(
remove_folded_structure = cls.remove_folded_structure,
recenter = cls.recenter,
)
#end def set_operations
def __init__(self,
axes = None,
scale = 1.,
elem = None,
pos = None,
elem_pos = None,
mag = None,
center = None,
kpoints = None,
kweights = None,
kgrid = None,
kshift = None,
permute = None,
units = None,
tiling = None,
rescale = True,
dim = 3,
operations = None,
background_charge = 0,
frozen = None,
bconds = None,
posu = None,
use_prim = None,
add_kpath = False,
symm_kgrid = False,
):
if isinstance(axes,str):
axes = np.array(axes.split(),dtype=float)
npe.reshape_inplace(axes, (dim, dim))
#end if
if center is None:
if axes is not None:
center = np.array(axes,dtype=float).sum(0)/2
else:
center = dim*[0]
#end if
#end if
if bconds is None or bconds=='periodic':
bconds = dim*['p']
#end if
if axes is None:
axes = []
bconds = []
#end if
if elem_pos is not None:
ep = np.array(elem_pos.split(),dtype=str)
npe.reshape_inplace(ep, (ep.size//(dim+1), (dim+1)))
elem = ep[:,0].ravel()
pos = ep[:,1:dim+1]
#end if
if elem is None:
elem = []
#end if
if posu is not None:
pos = posu
#end if
if pos is None:
pos = np.empty((0,dim))
#end if
if kshift is None:
kshift = 0,0,0
#end if
self.scale = 1.
self.units = units
self.dim = dim
self.center = np.array(center,dtype=float)
self.axes = np.array(axes,dtype=float)
self.set_bconds(bconds)
self.set_elem(elem)
self.set_pos(pos)
self.set_mag(mag)
self.set_frozen(frozen)
self.kpoints = np.empty((0,dim))
self.kweights = np.empty((0,))
self.background_charge = background_charge
self.remove_folded_structure()
if len(axes)==0:
self.kaxes=np.array([])
else:
self.kaxes=2*pi*inv(self.axes).T
#end if
if posu is not None:
self.pos_to_cartesian()
#end if
if use_prim is not None and use_prim is not False:
self.become_primitive(source=use_prim,add_kpath=add_kpath)
#end if
if tiling is not None:
self.tile(tiling,in_place=True)
#end if
if kpoints is not None:
self.add_kpoints(kpoints,kweights)
#end if
if kgrid is not None:
if not symm_kgrid:
self.add_kmesh(kgrid,kshift)
else:
self.add_symmetrized_kmesh(kgrid,kshift)
#end if
#end if
if rescale:
self.rescale(scale)
else:
self.scale = scale
#end if
if permute is not None:
self.permute(permute)
#end if
if operations is not None:
self.operate(operations)
#end if
#end def __init__
[docs]
def check_consistent(self,tol=1e-8,exit=True,message=False):
msg = ''
if self.has_axes():
kaxes = 2*pi*inv(self.axes).T
abs_diff = np.abs(self.kaxes-kaxes).sum()
if abs_diff>tol:
msg += (
"Direct and reciprocal space axes are not consistent.\n"
"axes present:\n"
"{0}\n"
"kaxes present:\n"
"{1}"
"\nConsistent kaxes:\n"
"{2}\n"
"Absolute difference: {3}\n".format(
self.axes, self.kaxes, kaxes,abs_diff
)
)
#end if
#end if
N = len(self.elem)
D = self.dim
pshape = (N,D)
if self.pos.shape!=pshape:
msg += (
"pos is not the right shape\n"
"pos shape: {}\n"
"Correct shape: {}\n".format(
self.pos.shape, pshape
)
)
#end if
if self.mag is not None and len(self.mag)!=N:
msg += (
"mag does not have the right length\n"
"mag length: {}\n"
"Correct length: {}\n".format(
self.mag, N
)
)
#end if
if self.frozen is not None and self.frozen.shape!=pshape:
msg += (
"frozen is not the right shape\n"
"frozen shape: {}\n"
"Correct shape: {}\n".format(
self.frozen.shape, pshape
)
)
#end if
consistent = len(msg)==0
if not consistent and exit:
self.error(msg)
#end if
if not message:
return consistent
else:
return consistent,msg
#end if
#end def check_consistent
[docs]
def set_axes(self,axes):
self.reset_axes(axes)
#end def set_axes
[docs]
def set_bconds(self,bconds):
self.bconds = np.array(tuple(bconds),dtype=str)
#end def bconds
[docs]
def set_elem(self,elem):
self.elem = np.array(elem,dtype=object)
#end def set_elem
[docs]
def set_pos(self,pos):
self.pos = np.array(pos,dtype=float)
if len(self.pos)!=len(self.elem):
self.error(
"Atomic positions must have same length as elem.\n"
"elem length: {}\n"
"Atomic positions length: {}\n".format(
len(self.elem), len(self.pos)
)
)
#end if
#end def set_pos
[docs]
def set_mag(self,mag=None):
if mag is None:
self.mag = None
else:
self.mag = np.array(mag,dtype=object)
if len(self.mag)!=len(self.elem):
self.error(
"Magnetic moments must have same length as elem.\n"
"elem length: {}\n"
"Magnetic moments length: {}\n".format(
len(self.elem), len(self.mag)
)
)
#end if
#end if
#end def set_mag
[docs]
def set_frozen(self,frozen=None):
if frozen is None:
self.frozen = None
else:
self.frozen = np.array(frozen,dtype=bool)
if self.frozen.shape!=self.pos.shape:
self.error(
"Frozen directions must have the same shape as positions.\n"
"Positions shape: {0}\n"
"Frozen directions shape: {1}".format(
self.pos.shape, self.frozen.shape
)
)
#end if
#end if
#end def set_frozen
[docs]
def size(self):
return len(self.elem)
#end def size
[docs]
def has_axes(self):
return len(self.axes)==self.dim
#end def has_axes
[docs]
def operate(self,operations):
for op in operations:
if op not in self.operations:
self.error(
"{0} is not a known operation\n"
"valid options are:\n"
" {1}".format(
op, list(self.operations.keys())
)
)
else:
self.operations[op](self)
#end if
#end for
#end def operate
[docs]
def has_tmatrix(self):
return 'tmatrix' in self and self.tmatrix is not None
#end def has_tmatrix
[docs]
def is_tiled(self):
return self.has_folded() and self.has_tmatrix()
#end def is_tiled
[docs]
def set_folded(self,folded):
self.set_folded_structure(folded)
#end def set_folded
[docs]
def remove_folded(self):
self.remove_folded_structure()
#end def remove_folded
[docs]
def has_folded(self):
return self.has_folded_structure()
#end def has_folded
[docs]
def set_folded_structure(self,folded):
if not isinstance(folded,Structure):
self.error(
"cannot set folded structure\n"
"folded structure must be an object with type Structure\n"
"received type: {0}".format(
folded.__class__.__name__
)
)
#end if
self.folded_structure = folded
if self.has_axes():
self.tmatrix = self.tilematrix(folded)
#end if
#end def set_folded_structure
[docs]
def remove_folded_structure(self):
self.folded_structure = None
if 'tmatrix' in self:
del self.tmatrix
#end if
#end def remove_folded_structure
[docs]
def has_folded_structure(self):
return self.folded_structure is not None
#end def has_folded_structure
[docs]
def group_atoms(self, folded=True) -> None:
"""Group the atoms by their element type, sorting in alphabetical order.
Parameters
----------
folded : bool, default=True
Optionally sort the folded structure, if it exists.
Examples
--------
>>> structure = Structure(
... elem = ["Be", "N", "C", "Be", "C", "O"],
... pos = np.array([
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... ], dtype=np.float64),
... )
>>> print(structure.elem)
['Be' 'N' 'C' 'Be' 'C' 'O']
>>> structure.group_atoms()
>>> print(structure.elem)
['Be' 'Be' 'C' 'C' 'N' 'O']
"""
if len(self.elem)>0:
order = self.elem.argsort()
if (self.elem!=self.elem[order]).any():
self.elem = self.elem[order]
self.pos = self.pos[order]
#end if
#end if
if self.folded_structure is not None and folded:
self.folded_structure.group_atoms(folded)
#end if
#end def group_atoms
[docs]
def rename(self, folded=True, **name_pairs) -> None:
"""Rename element names in a structure.
Parameters
----------
folded : bool, default=True
Rename elements in folded structure as well as tiled structure
(if there is no folded structure then this does nothing)
**name_pairs : dict[str]
A dictionary containing key:value pairs where the key is the old
element name and the value is the new element name.
Examples
--------
>>> structure = Structure(
... elem = ["N", "C", "O", "H"],
... pos = np.array([
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... [0, 0, 0],
... ], dtype=np.float64),
... )
>>> print(structure.elem)
['N' 'C' 'O' 'H']
>>> structure.rename(N="Dy", C="Er")
>>> print(structure.elem)
['Dy' 'Er' 'O' 'H']
"""
elem = self.elem
for old,new in name_pairs.items():
for i in range(len(self.elem)):
if old==elem[i]:
elem[i] = new
#end if
#end for
#end for
if self.folded_structure is not None and folded:
self.folded_structure.rename(folded=folded,**name_pairs)
#end if
#end def rename
[docs]
def reset_axes(self, axes: npt.ArrayLike | None = None) -> None:
"""Reset the structure's axes, k-space axes, and center.
Notes
-----
If ``axes`` is given, this function will remove the previous folded
structure (as new axes could invalidate the tiling) and then set
``self.axes`` to the new axes, update the k-space axes, and set the
center of the cell to be the (0.5, 0.5, 0.5) point.
If ``axes`` are not given (e.g. ``axes=None``, the default), then this
function will reuse the same axes as the current structure, reset the
k-space axes, and set the cell center at the (0.5, 0.5, 0.5) point.
Examples
--------
Providing axes will also make sure ``kaxes`` and ``center`` are consistent.
>>> structure = Structure(
... axes = np.array([
... [6.0, 0.0, 0.0],
... [0.0, 6.0, 0.0],
... [0.0, 0.0, 6.0],
... ]),
... )
>>> print(structure.axes)
[[6. 0. 0.]
[0. 6. 0.]
[0. 0. 6.]]
>>> print(structure.kaxes)
[[1.04719755 0. 0. ]
[0. 1.04719755 0. ]
[0. 0. 1.04719755]]
>>> print(structure.center)
[3. 3. 3.]
>>> structure.reset_axes(np.array([
... [12.0, 0.0, 0.0],
... [ 0.0, 12.0, 0.0],
... [ 0.0, 0.0, 12.0],
... ]))
>>> print(structure.axes)
[[12. 0. 0.]
[ 0. 12. 0.]
[ 0. 0. 12.]]
>>> print(structure.kaxes)
[[0.52359878 0. 0. ]
[0. 0.52359878 0. ]
[0. 0. 0.52359878]]
>>> print(structure.center)
[6. 6. 6.]
Alternatively, if you don't provide ``axes`` you can ensure that the
``kaxes`` and ``center`` of the structure are correct. Here we set
``center`` to be at the origin, then use ``reset_axes`` to correct it.
>>> structure = Structure(
... axes = np.array([
... [6.0, 0.0, 0.0],
... [0.0, 6.0, 0.0],
... [0.0, 0.0, 6.0],
... ]),
... center = np.array([0.0, 0.0, 0.0]),
... )
>>> print(structure.center)
[0. 0. 0.]
>>> structure.reset_axes()
>>> print(structure.center)
[3. 3. 3.]
"""
if axes is None:
axes = self.axes
else:
axes = np.array(axes)
self.remove_folded_structure()
#end if
self.axes = axes
self.kaxes = 2*pi*inv(axes).T
self.center = axes.sum(0)/2
#end def reset_axes
# test needed
[docs]
def adjust_axes(self,axes):
self.skew(dot(inv(self.axes),axes))
#end def adjust_axes
# test needed
[docs]
def reshape_axes(self,reshaping):
R = np.array(reshaping)
if np.abs(np.abs(det(R))-1)<1e-6:
self.axes = dot(self.axes,R)
else:
R = dot(inv(self.axes),R)
if np.abs(np.abs(det(R))-1)<1e-6:
self.axes = dot(self.axes,R)
else:
self.error(
"reshaping matrix must not change the volume\n"
" reshaping matrix:\n"
" {0}\n"
" volume change ratio: {1}".format(
R, np.abs(det(R))
)
)
#end if
#end if
#end def reshape_axes
[docs]
def write_axes(self) -> str:
"""Write the unit cell axes as a string.
Only implemented for ``self.dim=3``.
"""
c = ''
for a in self.axes:
c+='{0:12.8f} {1:12.8f} {2:12.8f}\n'.format(a[0],a[1],a[2])
#end for
return c
#end def write_axes
[docs]
def corners(self) -> npt.NDArray[np.float64]:
"""Calculate vectors corresponding to the 8 corners of the unit cell.
Only implemented for ``self.dim=3``.
"""
a = self.axes
c = np.array([(0,0,0),
a[0],
a[1],
a[2],
a[0]+a[1],
a[1]+a[2],
a[2]+a[0],
a[0]+a[1]+a[2],
])
return c
#end def corners
# test needed
[docs]
def miller_direction(self,h,k,l,normalize=False):
d = dot((h,k,l),self.axes)
if normalize:
d/=norm(d)
#end if
return d
#end def miller_direction
# test needed
[docs]
def miller_normal(self,h,k,l,normalize=False):
d = dot((h,k,l),self.kaxes)
if normalize:
d/=norm(d)
#end if
return d
#end def miller_normal
# test needed
[docs]
def project_plane(self,a1,a2,points=None):
# a1/a2: in plane vectors
if points is None:
points = self.pos
#end if
a1n = norm(a1)
a2n = norm(a2)
a1/=a1n
a2/=a2n
n = cross(a1,a2)
plane_coords = []
for p in points:
p -= dot(n,p)*n # project point into plane
c1 = dot(a1,p)/a1n
c2 = dot(a2,p)/a2n
plane_coords.append((c1,c2))
#end for
return np.array(plane_coords,dtype=float)
#end def project_plane
[docs]
def bounding_box(self,scale=1.0,minsize=None,mindist=0,box='tight',recenter=False):
pmin = self.pos.min(0)-mindist
pmax = self.pos.max(0)+mindist
pcenter = (pmax+pmin)/2
prange = pmax-pmin
if minsize is not None:
for i,pr in enumerate(prange):
prange[i] = max(minsize,prange[i])
#end for
#end if
if box=='tight':
axes = np.diag(prange)
elif box=='cubic' or box=='cube':
prmax = prange.max()
axes = np.diag((prmax,prmax,prmax))
elif isinstance(box,np.ndarray) or isinstance(box,list):
box = np.array(box)
if box.shape!=(3,3):
self.error(
"requested box must be 3-dimensional (3x3 axes)\n"
" you provided: "+str(box)+"\n"
" shape: "+str(box.shape)
)
#end if
binv = inv(box)
pu = dot(self.pos,binv)
pmin = pu.min(0)
pmax = pu.max(0)
pcenter = (pmax+pmin)/2
prange = pmax-pmin
axes = dot(np.diag(prange),box)
else:
self.error(
"invalid request for box\n"
" valid options are 'tight', 'cubic', or axes array (3x3)\n"
" you provided: "+str(box)
)
#end if
self.reset_axes(scale*axes)
self.slide(self.center-pcenter,recenter)
#end def bounding_box
[docs]
def center_molecule(self, rmin_edge = 1e-8):
"""Center a molecule in a unit cell, ensuring equal padding on all sides.
Parameters
----------
rmin_edge : float-like or array-like, default=1e-8
Minimum acceptable distance from the cell edges to any atom in the molecule, potentially specified per-axis.
Notes
-----
The default rmin_edge of 1e-8 is designed to just barely contain the molecule. It is
important to ensure that for larger, non-periodic molecular systems or atoms that the unit
cell size is sufficiently large to contain the charge density otherwise the molecular system
will interact with the neighboring periodic replica, resulting in unphysical electrostatic interactions
For reference, 5 Bohr contains >99.3% of the charge density around a Hydrogen atom.
"""
if hasattr(rmin_edge, '__len__'):
rmin_edge = np.abs(np.asarray(rmin_edge, dtype=float).flatten())
if rmin_edge.shape != (self.dim,):
self.error('rmin_edge must have length '+str(self.dim))
else:
rmin_edge = abs(float(rmin_edge))
#cell origin
corner = self.center - self.axes.sum(0) / 2
# Convert min distance to cell units
L = norm(self.axes, axis=0)
umin = rmin_edge / L
# Obtain current unit coordinates of the atoms
upos = np.dot(self.pos - corner, inv(self.axes))
# Align the molecule to contact the cell facets
upos -= upos.min(0)
# Align the molecule to be halfway between each cell facet
upos += (1-upos.max(0)) / 2
# Check that the molecule fits inside the box
utol = 1e-8 / L # numerical tolerance, 10 nano Bohr/Angstrom
upos_min = upos.min(0)
upos_max = upos.max(0)
lower_out = (upos_min < utol).any()
upper_out = (upos_max > (1.-utol)).any()
if lower_out or upper_out:
self.error('molecule does not fit in the cell')
# Check that the molecule is far enough from the cell edges
lower_out = (upos_min < umin).any()
upper_out = (upos_max > (1.-umin)).any()
if lower_out or upper_out:
self.error('molecule is too close to the cell edges')
# Transform back into real space coordinates
self.pos = np.dot(upos, self.axes) + corner
#end def center_molecule
# test needed
[docs]
def center_solid(self):
u = self.pos_unit()
du = (1-u.min(0)-u.max(0))/2
self.slide(dot(du,self.axes),recenter=False)
#end def center_solid
# test needed
[docs]
def permute(self,permutation):
dim = self.dim
P = np.empty((dim,dim),dtype=int)
if len(permutation)!=dim:
self.error(
" permutation vector must have {0} elements\n"
" you provided {1}".format(
dim, permutation
)
)
#end if
for i in range(dim):
p = permutation[i]
pv = np.zeros((dim,),dtype=int)
if p=='x' or p=='0':
pv[0] = 1
elif p=='y' or p=='1':
pv[1] = 1
elif p=='z' or p=='2':
pv[2] = 1
#end if
P[:,i] = pv[:]
#end for
self.center = dot(self.center,P)
if self.has_axes():
self.axes = dot(self.axes,P)
#end if
if len(self.pos)>0:
self.pos = dot(self.pos,P)
#end if
if len(self.kaxes)>0:
self.kaxes = dot(self.kaxes,P)
#end if
if len(self.kpoints)>0:
self.kpoints = dot(self.kpoints,P)
#end if
if self.folded_structure is not None:
self.folded_structure.permute(permutation)
#end if
#end def permute
# test needed
[docs]
def rotate_plane(self,plane,angle,units='degrees'):
self.pos = rotate_plane(plane,angle,self.pos,units)
if self.has_axes():
axes = rotate_plane(plane,angle,self.axes,units)
self.reset_axes(axes)
#end if
#end def rotate_plane
# test needed
[docs]
def upcast(self,DerivedStructure):
if not issubclass(DerivedStructure,Structure):
self.error(DerivedStructure.__name__,'is not derived from Structure')
#end if
ds = DerivedStructure()
for name,value in self.items():
ds[name] = deepcopy(value)
#end for
return ds
#end def upcast
# test needed
[docs]
def incorporate(self,other):
self.set_elem(list(self.elem)+list(other.elem))
self.pos=np.array(list(self.pos)+list(other.pos))
#end def incorporate
# test needed
[docs]
def clone_from(self,other):
if not isinstance(other,Structure):
self.error(
"cloning failed\n"
"can only clone from other Structure objects\n"
"received object of type: {0}".format(
other.__class__.__name__
)
)
#end if
o = other.copy()
self.__dict__ = o.__dict__
#end def clone_from
# test needed
[docs]
def add_atoms(self,elem,pos):
self.set_elem(list(self.elem)+list(elem))
self.pos=np.array(list(self.pos)+list(pos))
#end def add_atoms
[docs]
def is_open(self):
return not self.any_periodic()
#end def is_open
[docs]
def is_periodic(self):
return self.any_periodic()
#end def is_periodic
[docs]
def any_periodic(self):
has_cell = self.has_axes()
pbc = False
for bc in self.bconds:
pbc |= bc=='p'
#end if
periodic = has_cell and pbc
return periodic
#end def any_periodic
[docs]
def all_periodic(self):
has_cell = self.has_axes()
pbc = True
for bc in self.bconds:
pbc &= bc=='p'
#end if
periodic = has_cell and pbc
return periodic
#end def all_periodic
# test needed
[docs]
def distances(self,pos1=None,pos2=None):
if isinstance(pos1,Structure):
pos1 = pos1.pos
#end if
if pos2 is None:
if pos1 is None:
return sqrt((self.pos**2).sum(1))
else:
pos2 = self.pos
#end if
#end if
if len(pos1)!=len(pos2):
self.error('positions arrays are not the same length')
#end if
return sqrt(((pos1-pos2)**2).sum(1))
#end def distances
[docs]
def count_kshells(self, kcut, tilevec=[12, 12, 12], nkdig=10):
# check tilevec input
for nt in tilevec:
if nt % 2 != 0:
msg = 'tilevec must contain even integers'
msg += ' so that kgrid can be zero centered.'
Structure.class_error(msg, 'count_kshells')
#end if
#end for
origin = np.array([[0, 0, 0]])
axes = self.axes
raxes = 2*np.pi*inv(axes).T
kvecs = self.tile_points(origin, raxes, tilevec)
kvecs -= np.dot(tilevec, raxes)/2 # center around 0
kmags = norm(kvecs, axis=-1)
# make sure tilevec is sufficient for kcut
klimit = 0.5*kmags.max()
if kcut > klimit:
msg = 'kcut %3.2f > klimit=%3.2f\n' % (kcut, klimit)
msg += ' please increase tilevec to be safe.\n'
Structure.class_error(msg, 'count_kshells')
#end if
sel = (0<kmags) & (kmags<kcut)
ukmags = np.unique(kmags[sel].round(nkdig))
return len(ukmags)
#end def count_kshells
[docs]
def volume(self):
if not self.has_axes():
return None
else:
return np.abs(det(self.axes))
#end if
#end def volume
[docs]
def rwigner(self,nc=5):
if self.dim!=3:
self.error('rwigner is currently only implemented for 3 dimensions')
#end if
rmin = 1e90
n=np.empty((1,3))
rng = tuple(range(-nc,nc+1))
for k in rng:
for j in rng:
for i in rng:
if i!=0 or j!=0 or k!=0:
n[:] = i,j,k
rmin = min(rmin,.5*norm(dot(n,self.axes)))
#end if
#end for
#end for
#end for
return rmin
#end def rwigner
[docs]
def rinscribe(self):
if self.dim!=3:
self.error('rinscribe is currently only implemented for 3 dimensions')
#end if
radius = 1e99
dim=3
axes=self.axes
for d in range(dim):
i = d
j = (d+1)%dim
rc = cross(axes[i,:],axes[j,:])
radius = min(radius,.5*np.abs(det(axes))/norm(rc))
#end for
return radius
#end def rinscribe
# test needed
[docs]
def rwigner_cube(self,*args,**kwargs):
cube = Structure()
a = self.volume()**(1./3)
cube.set_axes([[a,0,0],[0,a,0],[0,0,a]])
return cube.rwigner(*args,**kwargs)
#end def rwigner_cube
# test needed
[docs]
def rinscribe_cube(self,*args,**kwargs):
cube = Structure()
a = self.volume()**(1./3)
cube.set_axes([[a,0,0],[0,a,0],[0,0,a]])
return cube.rinscribe(*args,**kwargs)
#end def rinscribe_cube
[docs]
def rmin(self):
return self.rwigner()
#end def rmin
[docs]
def rcell(self):
return self.rinscribe()
#end def rcell
# test needed
# scale invariant measure of deviation from cube shape
# based on deviation of face diagonals from cube
[docs]
def cube_deviation(self):
a = self.axes
dc = self.volume()**(1./3)*sqrt(2.)
d1 = np.abs(norm(a[0]+a[1])-dc)
d2 = np.abs(norm(a[1]+a[2])-dc)
d3 = np.abs(norm(a[2]+a[0])-dc)
d4 = np.abs(norm(a[0]-a[1])-dc)
d5 = np.abs(norm(a[1]-a[2])-dc)
d6 = np.abs(norm(a[2]-a[0])-dc)
return (d1+d2+d3+d4+d5+d6)/(6*dc)
#end def cube_deviation
# test needed
# apply volume preserving shear-removing transformations to cell axes
# resulting unsheared cell has orthogonal axes
# while remaining periodically correct
# note that the unshearing procedure is not unique
# it depends on the order of unshearing operations
[docs]
def unsheared_axes(self,axes=None,distances=False):
if self.dim!=3:
self.error('unsheared_axes is currently only implemented for 3 dimensions')
#end if
if axes is None:
axes = self.axes
#end if
dim=3
axbar = np.identity(dim)
axnew = np.array(axes,dtype=float)
dists = np.empty((dim,))
for d in range(dim):
d2 = (d+1)%dim
d3 = (d+2)%dim
n = cross(axnew[d2],axnew[d3]) #vector normal to 2 cell faces
axdist = dot(n,axes[d])/dot(n,axbar[d])
axnew[d] = axdist*axbar[d]
dists[d] = axdist
#end for
if not distances:
return axnew
else:
return axnew,dists
#end if
#end def unsheared_axes
# test needed
# vectors parallel to cell faces
# length of vectors is distance between parallel face planes
# note that the product of distances is not the cell volume in general
# see "unsheared_axes" function
# (e.g. a volume preserving shear may bring two face planes arbitrarily close)
[docs]
def face_vectors(self,axes=None,distances=False):
if axes is None:
axes = self.axes
#end if
fv = inv(axes).T
for d in range(len(fv)):
fv[d] /= norm(fv[d]) # face normals
#end for
dv = dot(axes,fv.T) # axis projections onto face normals
fv = dot(dv,fv) # face normals lengthened by plane separation
if not distances:
return fv
else:
return fv,np.diag(dv)
#end if
#end def face_vectors
# test needed
[docs]
def face_distances(self):
return self.face_vectors(distances=True)[1]
#end def face_distances
# test needed
[docs]
def rescale(self,scale):
self.scale *= scale
self.axes *= scale
self.pos *= scale
self.center *= scale
self.kaxes /= scale
self.kpoints/= scale
if self.folded_structure is not None:
self.folded_structure.rescale(scale)
#end if
#end def rescale
# test needed
[docs]
def stretch(self,s1,s2,s3):
if self.dim!=3:
self.error('stretch is currently only implemented for 3 dimensions')
#end if
d = np.diag((s1,s2,s3))
self.skew(d)
#end def stretch
# test needed
[docs]
def rotate(self,r,rp=None,passive=False,units="radians",check=True):
"""Arbitrary rotation of the structure.
Parameters
----------
r : array_like of float with shape (3,3) or array_like of float shape (3,) or str
If a 3x3 matrix, then code executes rotation consistent with this
matrix -- it is assumed that the matrix acts on a column-major
vector (eg, v'=Rv). If a three-dimensional array, then the operation
of the function depends on the input type of ``rp`` in the
following ways:
1. If ``rp`` is a scalar, then ``rp`` is assumed to be an angle and
a rotation of ``rp`` is made about the axis defined by ``r``.
2. If ``rp`` is a vector, then ``rp`` is assumed to be an axis and a
rotation is made such that ``r`` aligns with ``rp``.
3. If ``rp`` is a ``str``, then the rotation is such that ``r``
aligns with the axis given by the str
``('x', 'y', 'z', 'a0', 'a1', or 'a2')``.
If a ``str`` then the axis, ``r``, is defined by the input label
(e.g. 'x', 'y', 'z', 'a1', 'a2', or 'a3') and the operation of the
function depends on the input type of ``rp`` in the same ways as
above.
rp : array_like of float with shape (3) or str, optional
If a 3-dimensional vector is given, then ``rp`` is assumed to be an
axis and a rotation is made such that the axis ``r`` is aligned with
``rp``.
If a str, then ``rp`` is assumed to be an angle and a
rotation about the axis defined by ``r`` is made by an angle ``rp``.
If a str is given, then ``rp`` is assumed to be an axis defined by
the given label (e.g. 'x', 'y', 'z', 'a1', 'a2', or 'a3') and a
rotation is made such that the axis ``r`` is aligned with ``rp``.
passive : bool, default=False
If ``True``, perform a passive rotation.
If ``False``, perform an active rotation.
units : {"radians", "rad", "degrees", "deg"}, default="radians"
Units of ``rp``, if ``rp`` is given as an angle (scalar).
check : bool, default=True
Perform a check to verify rotation matrix is orthogonal.
"""
if rp is not None:
dirmap = dict(x=[1,0,0],y=[0,1,0],z=[0,0,1])
if isinstance(r,str):
if r[0]=='a': # r= 'a0', 'a1', or 'a2'
r = self.axes[int(r[1])]
else: # r= 'x', 'y', or 'z'
r = dirmap[r]
#end if
else:
r = np.array(r,dtype=float)
if len(r.shape)>1:
self.error('r must be given as a 1-d vector or string, if rp is not None')
#end if
#end if
if isinstance(rp,(int,float)):
if units=="radians" or units=="rad":
theta = float(rp)
else:
theta = float(rp)*np.pi/180.0
#end if
c = np.cos(theta)
s = np.sin(theta)
else:
if isinstance(rp,str):
if rp[0]=='a': # rp= 'a0', 'a1', or 'a2'
rp = self.axes[int(rp[1])]
else: # rp= 'x', 'y', or 'z'
rp = dirmap[rp]
#end if
else:
rp = np.array(rp,dtype=float)
#end if
# go from r,rp to r,theta
c = np.dot(r,rp)/norm(r)/norm(rp)
if np.abs(c-1)<1e-6:
s = 0.0
r = np.array([1,0,0])
else:
s = np.dot(cross(r,rp),cross(r,rp))/norm(r)/norm(rp)/norm(cross(r,rp))
r = cross(r,rp)/norm(cross(r,rp))
#end if
#end if
# make R from r,theta
R = [[ c+r[0]**2.0*(1.0-c), r[0]*r[1]*(1.0-c)-r[2]*s, r[0]*r[2]*(1.0-c)+r[1]*s],
[r[1]*r[0]*(1.0-c)+r[2]*s, c+r[1]**2.0*(1.0-c), r[1]*r[2]*(1.0-c)-r[0]*s],
[r[2]*r[0]*(1.0-c)-r[1]*s, r[2]*r[1]*(1.0-c)+r[0]*s, c+r[2]**2.0*(1.0-c)]]
else:
R = r
#end if
R = np.array(R,dtype=float)
if passive:
R = R.T
#end if
if check:
if not np.allclose(dot(R,R.T),np.identity(len(R))):
self.error('the function, rotate, must be given an orthogonal matrix')
#end if
#end if
self.matrix_transform(R)
#end def rotate
# test needed
#end if
#end def matrix_transform
# test needed
[docs]
def skew(self,skew):
"""
Arbitrary transformation matrix (row-major).
Parameters
----------
skew : `array_like, float, shape (3,3)`
Transform the structure using the matrix skew. It is assumed that
skew is in row-major form, i.e., it transforms a vector v as
v' = vT
"""
self.matrix_transform(skew.T)
#end def skew
# test needed
[docs]
def change_units(self,units,folded=True):
if units!=self.units:
scale = convert(1,self.units,units)
self.scale *= scale
self.axes *= scale
self.pos *= scale
self.center *= scale
self.kaxes /= scale
self.kpoints/= scale
self.units = units
#end if
if self.folded_structure is not None and folded:
self.folded_structure.change_units(units,folded=folded)
#end if
#end def change_units
# test needed
# insert sep space at loc along axis
# if sep<0, space is removed instead
[docs]
def cleave(self,axis,loc,sep=None,remove=False,tol=1e-6):
self.remove_folded_structure()
if isinstance(axis,int):
if sep is None:
self.error('separation induced by cleave must be provided')
#end if
v = self.face_vectors()[axis]
if isinstance(loc,float):
c = loc*v/norm(v)
#end if
else:
v = axis
#end if
c = np.array(c) # point on cleave plane
v = np.array(v) # normal vector to cleave plane, norm is cleave separation
if sep is not None:
v = np.abs(sep)*v/norm(v)
#end if
if norm(v)<tol:
return
#end if
vn = np.array(v/norm(v))
if sep is not None and sep<0:
v = -v # preserve the normal direction for atom identification, but reverse the shift direction
#end if
self.recorner() # want box contents to be static
if self.has_axes():
components = 0
dim = self.dim
axes = self.axes
for i in range(dim):
i2 = (i+1)%dim
i3 = (i+2)%dim
a2 = axes[i2]/norm(axes[i2])
a3 = axes[i3]/norm(axes[i3])
comp = np.abs(dot(a2,vn))+np.abs(dot(a3,vn))
if comp < 1e-6:
components+=1
iaxis = i
#end if
#end for
commensurate = components==1
if not commensurate:
self.error(
"cannot insert vacuum because cleave is incommensurate with the cell\n"
" cleave plane must be parallel to a cell face"
)
#end if
a = self.axes[iaxis]
#self.axes[iaxis] = (1.+dot(v,a)/dot(a,a))*a
self.axes[iaxis] = (1.+dot(v,v)/dot(v,a))*a
#end if
indices = []
pos = self.pos
for i in range(len(pos)):
p = pos[i]
comp = dot(p-c,vn)
if comp>0 or np.abs(comp)<tol:
pos[i] += v
indices.append(i)
#end if
#end for
if remove:
self.remove(indices)
#end if
#end def cleave
# test needed
[docs]
def translate(self,v):
v = np.array(v)
pos = self.pos
for i in range(len(pos)):
pos[i]+=v
#end for
self.center+=v
if self.folded_structure is not None:
self.folded_structure.translate(v)
#end if
#end def translate
# test needed
[docs]
def slide(self,v,recenter=True):
v = np.array(v)
pos = self.pos
for i in range(len(pos)):
pos[i]+=v
#end for
if recenter:
self.recenter()
#end if
if self.folded_structure is not None:
self.folded_structure.slide(v,recenter)
#end if
#end def slide
# test needed
[docs]
def zero_corner(self):
corner = self.center-self.axes.sum(0)/2
self.translate(-corner)
#end def zero_corner
# test needed
[docs]
def locate_simple(self,pos):
pos = np.array(pos)
if pos.shape==(self.dim,):
pos = [pos]
#end if
nn = nearest_neighbors(1,self.pos,pos)
return nn.ravel()
#end def locate_simple
# test needed
[docs]
def locate(self,identifiers,radii=None,exterior=False):
indices = None
if isinstance(identifiers,Structure):
cell = identifiers
indices = cell.inside(self.pos)
elif isinstance(identifiers,np.ndarray) and identifiers.dtype==bool:
indices = np.arange(len(self.pos))[identifiers]
elif isinstance(identifiers,int):
indices = [identifiers]
elif len(identifiers)>0 and isinstance(identifiers[0],int):
indices = identifiers
elif isinstance(identifiers,str):
atom = identifiers
indices = []
for i in range(len(self.elem)):
if self.elem[i]==atom:
indices.append(i)
#end if
#end for
elif len(identifiers)>0 and isinstance(identifiers[0],str):
indices = []
for atom in identifiers:
for i in range(len(self.elem)):
if self.elem[i]==atom:
indices.append(i)
#end if
#end for
#end for
#end if
if radii is not None or indices is None:
if indices is None:
pos = identifiers
else:
pos = self.pos[indices]
#end if
if isinstance(radii,float) or isinstance(radii,int):
radii = len(pos)*[radii]
elif radii is not None and len(radii)!=len(pos):
self.error(
"lengths of input radii and positions do not match\n"
" len(radii)={0}\n"
" len(pos)={1}".format(
len(radii), len(pos)
)
)
#end if
dtable = self.min_image_distances(pos)
indices = []
if radii is None:
for i in range(len(pos)):
indices.append(dtable[i].argmin())
#end for
else:
ipos = np.arange(len(self.pos))
for i in range(len(pos)):
indices.extend(ipos[dtable[i]<radii[i]])
#end for
#end if
#end if
if exterior:
indices = list(set(range(len(self.pos)))-set(indices))
#end if
return indices
#end def locate
[docs]
def freeze(self,identifiers=None,radii=None,exterior=False,negate=False,directions='xyz'):
if isinstance(identifiers,np.ndarray) and identifiers.shape==self.pos.shape and identifiers.dtype==bool:
if negate:
self.frozen = ~identifiers
else:
self.frozen = identifiers.copy()
#end if
return
#end if
if identifiers is None:
indices = np.arange(len(self.pos),dtype=int)
else:
indices = self.locate(identifiers,radii,exterior)
#end if
if len(indices)==0:
self.error('failed to select any atoms to freeze')
#end if
if isinstance(directions,str):
d = np.empty((3,),dtype=bool)
d[0] = 'x' in directions
d[1] = 'y' in directions
d[2] = 'z' in directions
directions = len(indices)*[d]
else:
directions = np.array(directions,dtype=bool)
#end if
if self.frozen is None:
self.frozen = np.zeros(self.pos.shape,dtype=bool)
#end if
frozen = self.frozen
i=0
if not negate:
for index in indices:
frozen[index] = directions[i]
i+=1
#end for
else:
for index in indices:
frozen[index] = directions[i]==False
i+=1
#end for
#end if
#end def freeze
[docs]
def is_frozen(self):
if self.frozen is None:
return np.zeros((len(self.pos),),dtype=bool)
else:
return self.frozen.sum(1)>0
#end if
#end def is_frozen
[docs]
def is_magnetic(self,tol=1e-8):
magnetic = False
if self.mag is not None:
for m in self.mag:
if m is not None and np.abs(m)>tol:
magnetic = True
break
#end if
#end for
#end if
return magnetic
#end def is_magnetic
# test needed
[docs]
def carve(self,identifiers):
indices = self.locate(identifiers)
if isinstance(identifiers,Structure):
sub = identifiers
sub.elem = self.elem[indices].copy()
sub.pos = self.pos[indices].copy()
else:
sub = self.copy()
sub.elem = self.elem[indices]
sub.pos = self.pos[indices]
#end if
sub.host_indices = np.array(indices)
return sub
#end def carve
# test needed
[docs]
def remove(self,identifiers):
indices = self.locate(identifiers)
keep = list(set(range(len(self.pos)))-set(indices))
erem = self.elem[indices]
prem = self.pos[indices]
self.elem = self.elem[keep]
self.pos = self.pos[keep]
self.remove_folded_structure()
return erem,prem
#end def remove
# test needed
[docs]
def replace(self,identifiers,elem=None,pos=None,radii=None,exterior=False):
indices = self.locate(identifiers,radii,exterior)
if isinstance(elem,Structure):
cell = elem
elem = cell.elem
pos = cell.pos
elif elem is None:
elem = self.elem
#end if
indices=np.array(indices)
elem=np.array(elem,dtype=object)
pos =np.array(pos)
nrem = len(indices)
nadd = len(pos)
if nadd<nrem:
ar = np.array(list(range(0,nadd)))
rr = np.array(list(range(nadd,nrem)))
self.elem[indices[ar]] = elem[:]
self.pos[indices[ar]] = pos[:]
self.remove(indices[rr])
elif nadd>nrem:
ar = np.array(list(range(0,nrem)))
er = np.array(list(range(nrem,nadd)))
self.elem[indices[ar]] = elem[ar]
self.pos[indices[ar]] = pos[ar]
ii = indices[ar[-1]]
self.set_elem( list(self.elem[0:ii])+list(elem[er])+list(self.elem[ii:]) )
self.pos = np.array( list(self.pos[0:ii])+list(pos[er])+list(self.pos[ii:]) )
else:
self.elem[indices] = elem[:]
self.pos[indices] = pos[:]
#end if
self.remove_folded_structure()
#end def replace
# test needed
[docs]
def replace_nearest(self,elem,pos=None):
if isinstance(elem,Structure):
cell = elem
elem = cell.elem
pos = cell.pos
#end if
nn = nearest_neighbors(1,self.pos,pos)
np = len(pos)
nps= len(self.pos)
d = np.empty((np,))
ip = np.array(list(range(np)))
ips= nn.ravel()
for i in ip:
j = ips[i]
d[i]=sqrt(((pos[i]-self.pos[j])**2).sum())
#end for
order = d.argsort()
ip = ip[order]
ips=ips[order]
replacable = np.empty((nps,))
replacable[:] = False
replacable[ips]=True
insert = []
last_replaced=nps-1
for n in range(np):
i = ip[n]
j = ips[n]
if replacable[j]:
self.pos[j] = pos[i]
self.elem[j]=elem[i]
replacable[j]=False
last_replaced = j
else:
insert.append(i)
#end if
#end for
insert=np.array(insert)
ii = last_replaced
if len(insert)>0:
self.set_elem( list(self.elem[0:ii])+list(elem[insert])+list(self.elem[ii:]) )
self.pos = np.array( list(self.pos[0:ii])+list(pos[insert])+list(self.pos[ii:]) )
#end if
self.remove_folded_structure()
#end def replace_nearest
# test needed
[docs]
def point_defect(self,identifiers=None,elem=None,dr=None):
if isinstance(elem,str):
elem = [elem]
if dr is not None:
dr = [dr]
#end if
#end if
if 'point_defects' not in self:
self.point_defects = obj()
#end if
point_defects = self.point_defects
ncenters = len(point_defects)
if identifiers is None:
index = ncenters
if index>=len(self.pos):
self.error(
"attempted to add a point defect at index {0}, which does not exist\n"
" for reference there are {1} atoms in the structure".format(
index, len(self.pos)
)
)
#end if
else:
indices = self.locate(identifiers)
if len(indices)>1:
self.error(
"{0} atoms were located by identifiers provided\n"
" a point defect replaces only a single atom\n"
" atom indices located: {1}".format(
len(indices), indices
)
)
#end if
index = indices[0]
#end if
if elem is None:
self.error(
"must supply substitutional elements comprising the point defect\n"
" expected a list or similar for input argument elem"
)
elif len(elem)>1 and dr is None:
self.error(
"must supply displacements (dr) since many atoms comprise the point defect"
)
elif dr is not None and len(elem)!=len(dr):
self.error(
"elem and dr must have the same length"
)
#end if
r = self.pos[index]
e = self.elem[index]
elem = np.array(elem)
pos = np.zeros((len(elem),len(r)))
if dr is None:
rc = r
for i in range(len(elem)):
pos[i] = r
#end for
else:
nrc = 0
rc = 0*r
dr = np.array(dr)
for i in range(len(elem)):
pos[i] = r + dr[i]
if norm(dr[i])>1e-5:
rc+=dr[i]
nrc+=1
#end if
#end for
if nrc==0:
rc = r
else:
rc = r + rc/nrc
#end if
#end if
point_defect = obj(
center = rc,
elem_replaced = e,
elem = elem,
pos = pos
)
point_defects.append(point_defect)
elist = list(self.elem)
plist = list(self.pos)
if len(elem)==0 or len(elem)==1 and elem[0]=='':
elist.pop(index)
plist.pop(index)
else:
elist[index] = elem[0]
plist[index] = pos[0]
for i in range(1,len(elem)):
elist.append(elem[i])
plist.append(pos[i])
#end for
#end if
self.set_elem(elist)
self.pos = np.array(plist)
self.remove_folded_structure()
#end def point_defect
# test needed
[docs]
def species(self,symbol=False):
if not symbol:
return set(self.elem)
else:
species_labels = set(self.elem)
species = set()
for e in species_labels:
is_elem, element = Elements.is_element(e, return_element=True)
species.add(element.symbol)
#end for
return species_labels,species
#end if
#end def species
# test needed
[docs]
def ordered_species(self,symbol=False):
speclab_set = set()
species_labels = []
if not symbol:
for e in self.elem:
if e not in speclab_set:
speclab_set.add(e)
species_labels.append(e)
#end if
#end for
return species_labels
else:
species = []
spec_set = set()
for e in self.elem:
is_elem, element = Elements.is_element(e, return_element=True)
symbol = element.symbol
if e not in speclab_set:
speclab_set.add(e)
species_labels.append(e)
#end if
if symbol not in spec_set:
spec_set.add(symbol)
species.append(symbol)
#end if
#end for
return species_labels,species
#end if
#end def ordered_species
# test needed
[docs]
def order_by_species(self,folded=False):
species = []
species_counts = []
elem_indices = []
spec_set = set()
for i in range(len(self.elem)):
e = self.elem[i]
if e not in spec_set:
spec_set.add(e)
species.append(e)
species_counts.append(0)
elem_indices.append([])
#end if
sindex = species.index(e)
species_counts[sindex] += 1
elem_indices[sindex].append(i)
#end for
elem_order = []
for elem_inds in elem_indices:
elem_order.extend(elem_inds)
#end for
self.reorder(elem_order)
if folded and self.folded_structure is not None:
self.folded_structure.order_by_species(folded)
#end if
return species,species_counts
#end def order_by_species
# test needed
[docs]
def reorder(self,order):
order = np.array(order)
self.elem = self.elem[order]
self.pos = self.pos[order]
#end def reorder
# test needed
# find layers parallel to a particular cell face
# layers are found by scanning a window of width dtol along the axis and counting
# the number of atoms within the window. window position w/ max number of atoms
# defines the layer. layer distance is the window position.
# the resolution of the scan is determined by dbin
# (axis length)/dbin is the number of fine bins
# dtol/dbin is the number of fine bins in the moving (boxcar) window
# plot=True: plot the layer histogram (fine hist and moving average)
# composition=True: return the composition of each layer (count of each species)
# returns an object containing indices of atoms in each layer by distance along axis
# example: structure w/ 3 layers of 4 atoms each at distances 3.0, 6.0, and 9.0 Angs.
# layers
# 3.0 = [ 0, 1, 2, 3 ]
# 6.0 = [ 4, 5, 6, 7 ]
# 9.0 = [ 8, 9,10,11 ]
[docs]
def layers(self,axis=0,dtol=0.03,dbin=0.01,plot=False,composition=False):
nbox = int(dtol/dbin)
if nbox%2==0:
nbox+=1
#end if
nwind = (nbox-1)//2
s = self.copy()
s.recenter()
vaxis = s.axes[axis]
daxis = norm(vaxis)
naxis = vaxis/daxis
dbin = dtol/nbox
nbins = int(np.ceil(daxis/dbin))
dbin = daxis/nbins
dbins = daxis*(np.arange(nbins)+.5)/nbins
dists = daxis*s.pos_unit()[:,axis]
hist = np.zeros((nbins,),dtype=int)
boxhist = np.zeros((nbins,),dtype=int)
ihist = obj()
iboxhist = obj()
index = 0
for d in dists:
ibin = int(np.floor(d/dbin))
hist[ibin]+=1
if ibin not in ihist:
ihist[ibin] = []
#end if
ihist[ibin].append(index)
index+=1
#end for
for ib in range(nbins):
for i in range(ib-nwind,ib+nwind+1):
n = hist[i%nbins]
if n>0:
boxhist[ib]+=n
if ib not in iboxhist:
iboxhist[ib] = []
#end if
iboxhist[ib].extend(ihist[i%nbins])
#end if
#end for
#end for
peaks = []
nlast=0
for ib in range(nbins):
n = boxhist[ib]
if nlast==0 and n>0:
pcur = []
peaks.append(pcur)
#end if
if n>0:
pcur.append(ib)
#end if
nlast = n
#end for
if boxhist[0]>0 and boxhist[-1]>0:
peaks[0].extend(peaks[-1])
peaks.pop()
#end if
layers = obj()
ip = []
for peak in peaks:
ib = peak[boxhist[peak].argmax()]
ip.append(ib)
pindices = iboxhist[ib]
ldist = dbins[ib] # distance is along an axis vector
faxis = self.face_vectors()[axis]
ldist = dot(ldist*naxis,faxis/norm(faxis))
layers[ldist] = np.array(pindices,dtype=int)
#end for
if plot:
plt.plot(dbins,boxhist,'b.-',label='boxcar histogram')
plt.plot(dbins,hist,'r.-',label='fine histogram')
plt.plot(dbins[ip],boxhist[ip],'rv',markersize=20)
plt.show()
plt.legend()
#end if
if not composition:
return layers
else:
return layers,self.layer_composition(layers)
#end if
#end def layers
# test needed
[docs]
def layer_composition(self,layers):
lcomp = obj()
for d,ind in layers.items():
comp = obj()
elem = self.elem[ind]
for e in elem:
if e not in comp:
comp[e] = 1
else:
comp[e] += 1
#end if
#end for
lcomp[d]=comp
#end for
return lcomp
#end def layer_composition
# test needed
[docs]
def shells(self,identifiers,radii=None,exterior=False,cumshells=False,distances=False,dtol=1e-6):
# get indices for 'core' and 'bulk'
# core is selected by identifiers, forms core for shells to be built around
# bulk is all atoms except for core
if identifiers=='point_defects':
if 'point_defects' not in self:
self.error('requested shells around point defects, but structure has no point defects')
#end if
core = []
for pd in self.point_defects:
core.append(pd.center)
#end for
core = np.array(core)
bulk_ind = self.locate(core,radii=dtol,exterior=True)
core_ind = self.locate(bulk_ind,exterior=True)
bulk = self.pos[bulk_ind]
else:
core_ind = self.locate(identifiers,radii,exterior)
bulk_ind = self.locate(core_ind,exterior=True)
core = self.pos[core_ind]
bulk = self.pos[bulk_ind]
#end if
bulk_ind = np.array(bulk_ind,dtype=int)
# build distance table between bulk and core
dtable = self.distance_table(bulk,core)
# find shortest distance for each bulk atom to any core atom and order by distance
dist = dtable.min(1)
ind = np.arange(len(bulk))
order = dist.argsort()
dist = dist[order]
ind = bulk_ind[ind[order]]
# find shells around the core
# the closest atom to the core starts the first shell and defines a shell distance
# other atoms are in the shell if within dtol distance of the first atom
# otherwise a new shell is started
ns = 0
ds = -1
shells = obj()
shells[ns] = list(core_ind) # first shell is all core atoms
dshells = [0.]
for n in range(len(dist)):
if np.abs(dist[n]-ds)>dtol:
shell = [ind[n]] # new shell starts with single atom
ns+=1
shells[ns] = shell
ds = dist[n] # shell distance is distance of this atom from core
dshells.append(ds)
else:
shell.append(ind[n])
#end if
#end for
dshells = np.array(dshells,dtype=float)
results = [shells]
if cumshells:
# assemble cumulative shells, ie cumshell[ns] = sum(shells[n],n=0 to ns)
cumshells = obj()
cumshells[0] = list(shells[0])
for ns in range(1,len(shells)):
cumshells[ns] = cumshells[ns-1]+shells[ns]
#end for
for ns,cshell in cumshells.items():
cumshells[ns] = np.array(cshell,dtype=int)
#end for
results.append(cumshells)
#end if
for ns,shell in shells.items():
shells[ns] = np.array(shell,dtype=int)
if distances:
results.append(dshells)
#end if
if len(results)==1:
results = results[0]
#end if
return results
#end def shells
# test needed
# find connected sets of atoms.
# indices is a list of atomic indices to consider (self.pos[indices] are their positions)
# atoms are considered connected if they are within rmax of each other
# order sets the maximum number of atoms in any connected graph
# order = 1 returns single atoms
# order = 2 returns dimers + order=1 results
# order = 3 returns trimers + order=2 results
# ...
# degree is explained w/ an example: a triangle of atoms 0,1,2 and a line of atoms 3,4,5 (3 & 5 are not neighbors)
# degree = False : returned object (cgraphs) has following structure:
# cgraphs[1] = [ (0,), (1,), (2,), (3,), (4,), (5,) ] # first order connected graphs (atoms)
# cgraphs[2] = [ (0,1), (0,2), (1,2), (3,4), (4,5) ] # second order connected graphs (dimers)
# cgraphs[3] = [ (0,1,2), (3,4,5) ] # third order connected graphs (trimers)
# degree = True : returned object (cgraphs) has following structure:
# cgraphs
# 1 # first order connected graphs (atoms)
# 0 # sum of vertex degrees is 0 (a single atom has no neighbors)
# (0,) = [ (0,), (1,), (2,), (3,), (4,), (5,) ] # graphs with vertex degree (0,)
# 2 # second order connected graphs (dimers)
# 2 # sum of vertex degrees is 2 (each atom is connected to 1 neighbor)
# (1,1) = [ (0,1), (0,2), (1,2), (3,4), (4,5) ] # graphs with vertex degree (1,1)
# 3 # third order connected graphs (trimers)
# 4 # sum of vertex degrees is 4 (2 atoms have 1 neighbor and 1 atom has 2)
# (1,1,2) = [ (3,5,4) ]
# 6 # sum of vertex degrees is 6 (each atom is connected to 2 others)
# (2,2,2) = [ (0,1,2) ] # graphs with vertex degree (2,2,2)
[docs]
def connected_graphs(self,order,indices=None,rmax=None,nmax=None,voronoi=False,degree=False,site_maps=False,**spec_max):
if indices is None:
indices = np.arange(len(self.pos),dtype=int)
pos = self.pos
else:
pos = self.pos[indices]
#end if
num_indices = len(indices)
neigh_table = []
actual_indices = None
if voronoi:
actual_indices = True
neighbors = self.voronoi_neighbors(indices,restrict=True,distance_ordered=False)
for nilist in neighbors:
neigh_table.append(nilist)
#end for
else:
actual_indices = False
elem = set(self.elem[indices])
spec = set(spec_max.keys())
if spec==elem or rmax is not None:
None
elif spec<elem and nmax is not None:
for e in elem:
if e not in spec:
spec_max[e] = nmax
#end if
#end for
#end if
# get neighbor table for subset of atoms specified by indices
nt,dt = self.neighbor_table(pos,pos,distances=True)
# determine how many neighbors to consider based on rmax (all are neighbors if rmax is None)
nneigh = np.zeros((num_indices,),dtype=int)
if len(spec_max)>0:
for n in range(num_indices):
nneigh[n] = min(spec_max[self.elem[n]],len(nt[n]))
#end for
elif rmax is None:
nneigh[:] = num_indices
else:
nneigh = (dt<rmax).sum(1)
#end if
for i in range(num_indices):
neigh_table.append(nt[i,1:nneigh[i]])
#end for
del nt,dt,nneigh,elem,spec,rmax
#end if
neigh_table = np.array(neigh_table,dtype=int)
# record which atoms are neighbors to each other
neigh_pairs = set()
if actual_indices:
for i in range(num_indices):
for ni in neigh_table[i]:
neigh_pairs.add((i,ni))
neigh_pairs.add((ni,i))
#end for
#end for
else:
for i in range(num_indices):
for ni in neigh_table[i]:
ii = indices[i]
jj = indices[ni]
neigh_pairs.add((ii,jj))
neigh_pairs.add((jj,ii))
#end for
#end for
#end if
# find the connected graphs
graphs_found = set() # map to contain tuples of connected atom's indices
cgraphs = obj()
for o in range(1,order+1): # organize by order
cgraphs[o] = []
#end for
if order>0:
cg = cgraphs[1]
for i in range(num_indices): # list of single atoms
gi = (i,) # graph indices
cg.append(gi) # add graph to graph list of order 1
graphs_found.add(gi) # add graph to set of all graphs
#end for
for o in range(2,order+1): # graphs of order o are found by adding all
cglast = cgraphs[o-1] # possible single neighbors to each graph of order o-1
cg = cgraphs[o]
for gilast in cglast: # all graphs of order o-1
for i in gilast: # all indices in each graph of order o-1
for ni in neigh_table[i]: # neighbors of selected atom in o-1 graph
gi = tuple(sorted(gilast+(ni,))) # new graph with neighbor added
if gi not in graphs_found and len(set(gi))==o: # add it if it is new and really is order o
graphs_found.add(gi) # add graph to set of all graphs
cg.append(gi) # add graph to graph list of order o
#end if
#end for
#end for
#end for
#end for
#end if
if actual_indices:
for o,cg in cgraphs.items():
cgraphs[o] = np.array(cg,dtype=int)
#end for
else:
# map indices back to actual atomic indices
for o,cg in cgraphs.items():
cgmap = []
for gi in cg:
#gi = np.array(gi)
gimap = tuple(sorted(indices[np.array(gi)]))
cgmap.append(gimap)
#end for
cgraphs[o] = np.array(sorted(cgmap),dtype=int)
#end for
#end if
# reorganize the graph listing by cluster and vertex degree, if desired
if degree:
#degree_map = obj()
cgraphs_deg = obj()
for o,cg in cgraphs.items():
dgo = obj()
cgraphs_deg[o] = dgo
for gi in cg:
di = np.zeros((o,),dtype=int)
for m in range(o):
i = gi[m]
for n in range(m+1,o):
j = gi[n]
if (i,j) in neigh_pairs:
di[m]+=1
di[n]+=1
#end if
#end for
#end for
d = int(di.sum())
dorder = di.argsort()
di = tuple(di[dorder])
gi = tuple(np.array(gi)[dorder])
if d not in dgo:
dgo[d]=obj()
#end if
dgd = dgo[d]
if di not in dgd:
dgd[di] = []
#end if
dgd[di].append(gi)
#degree_map[gi] = d,di
#end for
for dgd in dgo:
for di,dgi in dgd.items():
dgd[di]=np.array(sorted(dgi),dtype=int)
#end for
#end for
#end for
cgraphs = cgraphs_deg
#end if
if not site_maps:
return cgraphs
else:
cmaps = obj()
if not degree:
for order,og in cgraphs.items():
cmap = obj()
for slist in og:
for s in slist:
if s not in cmap:
cmap[s] = obj()
#end if
cmap[s].append(slist)
#end for
#end for
cmaps[order] = cmap
#end for
else:
for order,og in cgraphs.items():
for total_degree,tg in og.items():
for local_degree,lg in tg.items():
cmap = obj()
for slist in lg:
n=0
for s in slist:
d = local_degree[n]
if s not in cmap:
cmap[s] = obj()
#end if
if d not in cmap[s]:
cmap[s][d] = obj()
#end if
cmap[s][d].append(slist)
n+=1
#end for
#end for
cmaps.add_attribute_path((order,total_degree,local_degree),cmap)
#end for
#end for
#end for
#end if
return cgraphs,cmaps
#end if
#end def connected_graphs
# test needed
# returns connected graphs that are rings up to the requested order
# rings are constructed by pairing lines that share endpoints
# all vertices of a ring have degree two
[docs]
def ring_graphs(self,order,**kwargs):
# get all half order connected graphs
line_order = order//2+order%2+1
cgraphs = self.connected_graphs(line_order,degree=True,site_maps=False,**kwargs)
# collect half order graphs that are lines
lgraphs = obj()
for o in range(2,line_order+1):
total_degree = 2*o-2
vertex_degree = tuple([1,1]+(o-2)*[2])
lg = None
if o in cgraphs:
cg = cgraphs[o]
if total_degree in cg:
dg = cg[total_degree]
if vertex_degree in dg:
lg = dg[vertex_degree]
#end if
#end if
#end if
if lg is not None:
lg_end = obj()
for gi in lg:
end_key = tuple(sorted(gi[0:2])) # end points
if end_key not in lg_end:
lg_end[end_key] = []
#end if
lg_end[end_key].append(tuple(gi))
#end for
lgraphs[o] = lg_end
#end if
#end for
# contruct rings from lines that share endpoints
rgraphs = obj()
for o in range(3,order+1):
o1 = o/2+1 # split half order for odd, same for even,
o2 = o1+o%2
lg1 = lgraphs.get_optional(o1,None) # sets of half order lines
lg2 = lgraphs.get_optional(o2,None)
if lg1 is not None and lg2 is not None:
rg = []
rset = set()
for end_key,llist1 in lg1.items(): # list of lines sharing endpoints
if end_key in lg2:
llist2 = lg2[end_key] # second list of lines sharing endpoints
for gi1 in llist1: # combine line pairs into rings
for gi2 in llist2:
ri = tuple(sorted(set(gi1+gi2[2:]))) # ring indices
if ri not in rset and len(ri)==o: # exclude repeated lines or rings
rg.append(ri)
rset.add(ri)
#end if
#end for
#end for
#end if
#end for
rgraphs[o] = np.array(sorted(rg),dtype=int)
#end if
#end for
return rgraphs
#end def ring_graphs
# test needed
# find the centroid of a set of points/atoms in min image convention
[docs]
def min_image_centroid(self,points=None,indices=None):
if indices is not None:
points = self.pos[indices]
elif points is None:
self.error('points or images must be provided to min_image_centroid')
#end if
p = np.array(points,dtype=float)
cprev = p[0]+1e99
c = p[0]
while(norm(c-cprev)>1e-8):
p = self.cell_image(p,center=c)
cprev = c
c = p.mean(axis=0)
#end def min_image_centroid
return c
#end def min_image_centroid
# test needed
# find min image centroids of multiple sets of points/atoms
[docs]
def min_image_centroids(self,points=None,indices=None):
cents = []
if points is not None:
for p in points:
cents.append(self.min_image_centroid(p))
#end for
elif indices is not None:
for ind in indices:
cents.append(self.min_image_centroid(indices=ind))
#end for
else:
self.error('points or images must be provided to min_image_centroid')
#end if
return np.array(cents,dtype=float)
#end def min_image_centroids
[docs]
def min_image_vectors(self,points=None,points2=None,axes=None,pairs=True):
if points is None:
points = self.pos
elif isinstance(points,Structure):
points = points.pos
#end if
if axes is None:
axes = self.axes
#end if
axinv = inv(axes)
points = np.array(points)
single = points.shape==(self.dim,)
if single:
points = [points]
#end if
if points2 is None:
points2 = self.pos
elif isinstance(points2,Structure):
points2 = points2.pos
elif points2.shape==(self.dim,):
points2 = [points2]
#end if
npoints = len(points)
npoints2 = len(points2)
if pairs:
vtable = np.empty((npoints,npoints2,self.dim),dtype=float)
i=-1
for p in points:
i+=1
j=-1
for pp in points2:
j+=1
u = dot(pp-p,axinv)
vtable[i,j] = dot(u-np.floor(u+.5),axes)
#end for
#end for
result = vtable
else:
if npoints!=npoints2:
self.error(
"cannot create one to one minimum image vectors, point sets differ in length\n"
" npoints1 = {0}\n"
" npoints2 = {1}".format(
npoints, npoints2
)
)
#end if
vectors = np.empty((npoints,self.dim),dtype=float)
n = 0
for p in points:
pp = points2[n]
u = dot(pp-p,axinv)
vectors[n] = dot(u-np.floor(u+.5),axes)
n+=1
#end for
result = vectors
#end if
return result
#end def min_image_vectors
[docs]
def min_image_distances(self,points=None,points2=None,axes=None,vectors=False,pairs=True):
vtable = self.min_image_vectors(points,points2,axes,pairs=pairs)
rdim = len(vtable.shape)-1
dtable = sqrt((vtable**2).sum(rdim))
if not vectors:
return dtable
else:
return dtable,vtable
#end if
#end def min_image_distances
[docs]
def distance_table(self,points=None,points2=None,axes=None,vectors=False):
return self.min_image_distances(points,points2,axes,vectors)
#end def distance_table
[docs]
def vector_table(self,points=None,points2=None,axes=None):
return self.min_image_vectors(points,points2,axes)
#end def vector_table
[docs]
def neighbor_table(self,points=None,points2=None,axes=None,distances=False,vectors=False):
dtable,vtable = self.min_image_distances(points,points2,axes,vectors=True)
ntable = np.empty(dtable.shape,dtype=int)
for i in range(len(dtable)):
ntable[i] = dtable[i].argsort()
#end for
results = [ntable]
if distances:
for i in range(len(dtable)):
dtable[i] = dtable[i][ntable[i]]
#end for
results.append(dtable)
#end if
if vectors:
for i in range(len(vtable)):
vtable[i] = vtable[i][ntable[i]]
#end for
results.append(vtable)
#end if
if len(results)==1:
results = results[0]
#end if
return results
#end def neighbor_table
# test needed
[docs]
def min_image_norms(self,points,norms):
if isinstance(norms,int) or isinstance(norms,float):
norms = [norms]
#end if
vtable = self.min_image_vectors(points)
rdim = len(vtable.shape)-1
nout = []
for p in norms:
nout.append( ((np.abs(vtable)**p).sum(rdim))**(1./p) )
#end for
if len(norms)==1:
nout = nout[0]
#end if
return nout
#end def min_image_norms
# test needed
# get all neighbors according to contacting voronoi polyhedra in PBC
[docs]
def voronoi_neighbors(self,indices=None,restrict=False,distance_ordered=True):
if indices is None:
indices = np.arange(len(self.pos))
#end if
indices = set(indices)
# make a new version of this (small cell)
sn = self.copy()
sn.recenter()
# tile a large cell periodically
d = 3
t = tuple(np.zeros((d,),dtype=int)+3)
ss = sn.tile(t)
ss.recenter(sn.center)
# get nearest neighbor index pairs in the large cell
neigh_pairs = voronoi_neighbors(ss.pos)
# create a mapping from large to small indices
large_to_small = 3**d*list(range(len(self.pos)))
# find the neighbor pairs in the small cell
neighbors = obj()
small_inds = set(ss.locate(sn.pos))
for n in range(len(neigh_pairs)):
i,j = neigh_pairs[n,:]
if i in small_inds or j in small_inds: # pairs w/ at least one in cell image
i = large_to_small[i] # mapping to small cell indices
j = large_to_small[j]
if not restrict or (i in indices and j in indices): # restrict to orig index set
if i not in neighbors:
neighbors[i] = [j]
else:
neighbors[i].append(j)
#ned if
if j not in neighbors:
neighbors[j] = [i]
else:
neighbors[j].append(i)
#end if
#end if
#end if
#end for
# remove any duplicates and order by distance
if distance_ordered:
dt = self.distance_table()
for i,ni in neighbors.items():
ni = np.array(list(set(ni)),dtype=int)
di = dt[i,ni]
order = di.argsort()
neighbors[i] = ni[order]
#end for
else: # just remove duplicates
for i,ni in neighbors.items():
neighbors[i] = np.array(list(set(ni)),dtype=int)
#end for
#end if
return neighbors
#end def voronoi_neighbors
[docs]
def voronoi_vectors(self,indices=None,restrict=None):
ni = self.voronoi_neighbors(indices,restrict)
vt = self.vector_table()
vv = obj()
for i,vi in ni.items():
vv[i] = vt[i,vi]
#end for
return vv
#end def voronoi_vectors
[docs]
def voronoi_distances(self,indices=None,restrict=False):
vv = self.voronoi_vectors(indices,restrict)
vd = obj()
for i,vvi in vv.items():
vd[i] = norm(vvi,axis=1)
#end for
return vd
#end def voronoi_distances
[docs]
def voronoi_radii(self,indices=None,restrict=None):
vd = self.voronoi_distances(indices,restrict)
vr = obj()
for i,vdi in vd.items():
vr[i] = vdi.min()/2
#end for
return vr
#end def voronoi_radii
[docs]
def voronoi_species_radii(self):
vr = self.voronoi_radii()
vsr = obj()
for i,r in vr.items():
e = self.elem[i]
if e not in vsr:
vsr[e] = r
else:
vsr[e] = min(vsr[e],r)
#end if
#end for
return vsr
#end def voronoi_species_radii
# test needed
# get nearest neighbors according to constraints (voronoi, max distance, coord. number)
[docs]
def nearest_neighbors(self,indices=None,rmax=None,nmax=None,restrict=False,voronoi=False,distances=False,**spec_max):
if indices is None:
indices = np.arange(len(self.pos))
#end if
elem = set(self.elem[indices])
spec = set(spec_max.keys())
if spec==elem or rmax is not None or voronoi:
None
elif spec<elem and nmax is not None:
for e in elem:
if e not in spec:
spec_max[e] = nmax
#end if
#end for
else:
self.error(
"must specify nmax for all species\n"
" species present: {0}\n"
" you only provided nmax for these species: {1}".format(
sorted(elem), sorted(spec)
)
)
#end if
pos = self.pos[indices]
if not restrict:
pos2 = self.pos
else:
pos2 = pos
#end if
if voronoi:
neighbors = self.voronoi_neighbors(indices=indices,restrict=restrict)
dt = self.distance_table(pos,pos2)[:,1:]
else:
nt,dt = self.neighbor_table(pos,pos2,distances=True)
dt=dt[:,1:]
nt=nt[:,1:]
neighbors = list(nt)
#end if
for i in range(len(indices)):
neighbors[i] = indices[neighbors[i]]
#end for
dist = list(dt)
if rmax is None:
for i in range(len(indices)):
nn = neighbors[i]
dn = dist[i]
e = self.elem[indices[i]]
if e in spec_max:
smax = spec_max[e]
if len(nn)>smax:
neighbors[i] = nn[:smax]
dist[i] = dn[:smax]
#end if
#end if
#end for
else:
for i in range(len(indices)):
neighbors[i] = neighbors[i][dt[i]<rmax]
#end for
#end if
if not distances:
return neighbors
else:
return neighbors,dist
#end if
#end def nearest_neighbors
# test needed
# determine local chemical coordination limited by constraints
[docs]
def chemical_coordination(self,indices=None,nmax=None,rmax=None,restrict=False,voronoi=False,neighbors=False,distances=False,**spec_max):
if indices is None:
indices = np.arange(len(self.pos))
#end if
if not distances:
neigh = self.nearest_neighbors(indices=indices,nmax=nmax,rmax=rmax,restrict=restrict,voronoi=voronoi,**spec_max)
else:
neigh,dist = self.nearest_neighbors(indices=indices,nmax=nmax,rmax=rmax,restrict=restrict,voronoi=voronoi,distances=True,**spec_max)
#end if
neigh_elem = []
for i in range(len(indices)):
neigh_elem.extend(self.elem[neigh[i]])
#end for
chem_key = tuple(sorted(set(neigh_elem)))
chem_coord = np.zeros((len(indices),len(chem_key)),dtype=int)
for i in range(len(indices)):
counts = np.zeros((len(chem_key),),dtype=int)
nn = list(self.elem[neigh[i]])
for n in range(len(counts)):
chem_coord[i,n] = nn.count(chem_key[n])
#end for
#end for
chem_map = obj()
i=0
for coord in chem_coord:
coord = tuple(coord)
if coord not in chem_map:
chem_map[coord] = [indices[i]]
else:
chem_map[coord].append(indices[i])
#end if
i+=1
#end for
for coord,ind in chem_map.items():
chem_map[coord] = np.array(ind,dtype=int)
#end for
results = [chem_key,chem_coord,chem_map]
if neighbors:
results.append(neigh)
#end if
if distances:
results.append(dist)
#end if
return results
#end def chemical_coordination
# test needed
[docs]
def rcore_max(self,units=None):
nt,dt = self.neighbor_table(self.pos,distances=True)
d = dt[:,1]
rcm = d.min()/2
if units is not None:
rcm = convert(rcm,self.units,units)
#end if
return rcm
#end def rcore_max
# test needed
[docs]
def cell_image(self,p,center=None):
pos = np.array(p,dtype=float)
if center is None:
c = self.center.copy()
else:
c = np.array(center,dtype=float)
#end if
axes = self.axes
axinv = inv(axes)
for i in range(len(pos)):
u = dot(pos[i]-c,axinv)
pos[i] = dot(u-np.floor(u+.5),axes)+c
#end for
return pos
#end def cell_image
# test needed
[docs]
def center_distances(self,points,center=None):
if center is None:
c = self.center.copy()
else:
c = np.array(center,dtype=float)
#end if
points = self.cell_image(points,center=c)
for i in range(len(points)):
points[i] -= c
#end for
return sqrt((points**2).sum(1))
#end def center_distances
# test needed
[docs]
def recenter(self, center=None):
"""Center atoms around a new provided center of the unit cell, or if a new
center is not provided then use the (0.5 0.5 0.5) point of the unit cell.
"""
if center is not None:
self.center = np.array(center, dtype=float)
#end if
self.pos = recenter_points(pos=self.pos, center=self.center, axes=self.axes)
self.recenter_k()
#end def recenter
# test needed
[docs]
def recenter_k(self,kpoints=None,kaxes=None,kcenter=None,remove_duplicates=False):
"""Center k-points around the provided center of k-space.
Parameters
----------
kpoints : NDArray | None, default = None
Array of *N* k-points to center with shape (*N*,3)
kaxes : NDArray | None, default = None
Array of axes describing the k-space cell with shape (3,3)
kcenter : NDArray | None, default = None
Array containing a custom center of k-space with shape (3,3)
remove_duplicates : bool, default = False
Remove k-points that have a distance less than 1e-8
"""
use_self = kpoints is None
if use_self:
kpoints=self.kpoints
#end if
if kaxes is None:
kaxes=self.kaxes
#end if
if len(kpoints)>0:
axes = kaxes
if kcenter is None:
center = axes.sum(0)/2
else:
center = np.array(kcenter)
kpoints = recenter_points(pos=kpoints, center=center, axes=axes)
if remove_duplicates:
inside = self.inside(kpoints,axes,center)
kpoints = kpoints[inside]
nkpoints = len(kpoints)
unique = np.empty((nkpoints,),dtype=bool)
unique[:] = True
nn = nearest_neighbors(1,kpoints)
if nkpoints>1:
npe.reshape_inplace(nn, (nkpoints,))
dist = self.distances(kpoints,kpoints[nn])
tol = 1e-8
duplicates = np.arange(nkpoints)[dist<tol]
for i in duplicates:
if unique[i]:
for j in duplicates:
if sqrt(((kpoints[i]-kpoints[j])**2).sum(1))<tol:
unique[j] = False
#end if
#end for
#end if
#end for
#end if
kpoints = kpoints[unique]
#end if
#end if
if use_self:
self.kpoints = kpoints
else:
return kpoints
#end if
#end def recenter_k
# test needed
[docs]
def recorner(self, center = None):
"""Center atoms around the origin of the cell.
Parameters
----------
center : NDArray, default = self.center
Position of the center of the cell.
Notes
-----
If the user supplies ``center``, then this will modify
``self.center`` to reflect that change.
"""
if center is not None:
self.center = np.array(center, dtype=float)
pos = self.pos
axes = self.axes
corner = self.center - axes.sum(0)/2
self.pos = recenter_points(pos=pos, center=corner, axes=axes)
#end def recorner
# test needed
[docs]
def inside(self,pos,axes=None,center=None,tol=1e-8,separate=False):
if axes is None:
axes=self.axes
#end if
if center is None:
center=self.center
#end if
axes = np.array(axes)
center = np.array(center)
inside = []
surface = []
su = []
axinv = inv(axes)
for i in range(len(pos)):
u = dot(pos[i]-center,axinv)
umax = np.abs(u).max()
if np.abs(umax-.5)<tol:
surface.append(i)
su.append(u)
elif umax<.5:
inside.append(i)
#end if
#end for
npos,dim = pos.shape
drange = list(range(dim))
n = len(surface)
i=0
while i<n:
j=i+1
while j<n:
du = np.abs(su[i]-su[j])
match = False
for d in drange:
match = match or np.abs(du[d]-1.)<tol
#end for
if match:
surface[j]=surface[-1]
surface.pop()
su[j]=su[-1]
su.pop()
n-=1
else:
j+=1
#end if
#end while
i+=1
#end while
if not separate:
inside+=surface
return inside
else:
return inside,surface
#end if
#end def inside
[docs]
def tile(self,*td,**kwargs):
in_place = kwargs.pop('in_place',False)
check = kwargs.pop('check',False)
dim = self.dim
if len(td)==1:
if isinstance(td[0],int):
tiling = dim*[td[0]]
else:
tiling = td[0]
#end if
else:
tiling = td
#end if
tiling = np.array(tiling)
matrix_tiling = tiling.shape == (dim,dim)
tilematrix,tilevector = reduce_tilematrix(tiling)
ncells = int(np.round( np.abs(det(tilematrix)) ))
if ncells==1 and np.abs(tilematrix-np.identity(self.dim)).sum()<1e-1:
if in_place:
return self
else:
return self.copy()
#end if
#end if
self.recenter()
elem = np.array(ncells*list(self.elem))
pos = self.tile_points(self.pos,self.axes,tilematrix,tilevector)
axes = dot(tilematrix,self.axes)
center = axes.sum(0)/2
kaxes = dot(inv(tilematrix.T),self.kaxes)
kpoints = np.array(self.kpoints)
kweights = np.array(self.kweights)
mag = None
frozen = None
if self.mag is not None:
mag = ncells*list(self.mag)
#end if
if self.frozen is not None:
frozen = ncells*list(self.frozen)
#end if
ts = self.copy()
ts.center = center
ts.set_elem(elem)
ts.axes = axes
ts.pos = pos
ts.mag = mag
ts.kaxes = kaxes
ts.kpoints = kpoints
ts.kweights = kweights
ts.set_mag(mag)
ts.set_frozen(frozen)
ts.background_charge = ncells*self.background_charge
ts.recenter()
ts.unique_kpoints()
if self.is_tiled():
ts.tmatrix = dot(tilematrix,self.tmatrix)
ts.folded_structure = self.folded_structure.copy()
else:
ts.tmatrix = tilematrix
ts.folded_structure = self.copy()
#end if
if in_place:
self.clear()
self.transfer_from(ts)
ts = self
#end if
if check:
ts.check_tiling()
#end if
return ts
#end def tile
[docs]
def tile_points(self,points,axes,tilemat,tilevec=None):
if tilevec is None:
tilemat,tilevec = reduce_tilematrix(tilemat)
#end if
if not isinstance(tilemat,np.ndarray):
tilemat = np.array(tilemat)
#end if
matrix_tiling = np.abs(tilemat-np.diag(np.diag(tilemat))).sum()>0.1
if not matrix_tiling:
return self.tile_points_simple(points,axes,np.diag(np.abs(tilemat)))
else:
if not isinstance(axes,np.ndarray):
axes = np.array(axes)
#end if
if not isinstance(tilevec,np.ndarray):
tilevec = np.array(tilevec)
#end if
dim = len(axes)
npoints = len(points)
ntpoints = npoints*int(np.round(np.abs(det(tilemat))))
if tilevec.size==dim:
npe.reshape_inplace(tilevec, (1, dim))
#end if
taxes = dot(tilemat,axes)
success = False
for tvec in tilevec:
tpoints = self.tile_points_simple(points,axes,tvec)
tpoints,weights,pmap = self.unique_points_fast(tpoints,taxes)
if len(tpoints)==ntpoints:
success = True
break
#end if
#end for
if not success:
tpoints = self.tile_points_brute(points,axes,tilemat)
tpoints,weights,pmap = self.unique_points_fast(tpoints,taxes)
if len(tpoints)!=ntpoints:
self.error('brute force tiling failed')
#end if
#end if
#end if
return tpoints
#end def tile_points
[docs]
def tile_points_simple(self,points,axes,tilevec):
if not isinstance(points,np.ndarray):
points = np.array(points)
#end if
if not isinstance(tilevec,np.ndarray):
tilevec = np.array(tilevec)
#end if
if not isinstance(axes,np.ndarray):
axes = np.array(axes)
#end if
if len(points.shape)==1:
npoints,dim = len(points),1
else:
npoints,dim = points.shape
#end if
t = tilevec
ti = np.array(np.around(t),dtype=int)
noninteger = np.abs(t-ti).sum()>1e-6
if noninteger:
tp = t.prod()
if np.abs(tp-int(tp))>1e-6:
self.error(
"tiling vector does not correspond to an integer volume change\n"
"tiling vector: {0}\n"
"volume change: {1} {2} {3}".format(
tilevec,tilevec.prod(),npoints,int(npoints)
)
)
#end if
t = np.array(np.ceil(t),dtype=int)+1
else:
t = ti
#end if
if t.min()<0:
self.error(
"tiling vector cannot be negative\n"
"tiling vector provided: {}".format(t)
)
#end if
ntpoints = npoints*int(np.round( t.prod() ))
if ntpoints==0:
tpoints = np.array([])
else:
tpoints = np.empty((ntpoints,dim))
ns=0
ne=npoints
for k in range(t[2]):
for j in range(t[1]):
for i in range(t[0]):
v = dot(np.array([[i,j,k]]),axes)
for d in range(dim):
tpoints[ns:ne,d] = points[:,d]+v[0,d]
#end for
ns+=npoints
ne+=npoints
#end for
#end for
#end for
#end if
return tpoints
#end def tile_points_simple
[docs]
def tile_points_brute(self,points,axes,tilemat):
tcorners = [[0,0,0],
[1,0,0],
[0,1,0],
[0,0,1],
[0,1,1],
[1,0,1],
[1,1,0],
[1,1,1]]
tcorners = dot(tcorners,tilemat)
tmin = tcorners.min(axis=0)
tmax = tcorners.max(axis=0)
tilevec = tmax-tmin
tpoints = self.tile_points_simple(points,axes,tilevec)
return tpoints
#end def tile_points_brute
[docs]
def opt_tilematrix(self,*args,**kwargs):
return optimal_tilematrix(self.axes,*args,**kwargs)
#end def opt_tilematrix
[docs]
def tile_opt(self,*args,**kwargs):
Topt,ropt = self.opt_tilematrix(*args,**kwargs)
return self.tile(Topt)
#end def tile_opt
[docs]
def check_tiling(self,tol=1e-6,exit=True):
msg = ''
if not self.is_tiled():
return msg
#end if
msgs = []
st = self
s = self.folded_structure
nt = len(st.pos)
n = len(s.pos)
if nt%n!=0:
msgs.append('tiled atom count does is not divisible by untiled atom count')
#end if
vratio = st.volume()/s.volume()
if np.abs(vratio-float(nt)/n)>tol:
msgs.append('tiled/untiled volume ratio does not match tiled/untiled atom count ratio')
#end if
if np.abs(vratio-np.abs(det(st.tmatrix)))>tol:
msgs.append('tiled/untiled volume ratio does not match tiling matrix determinant')
#end if
p,w,pmap = self.unique_points_fast(st.pos,st.axes)
if len(p)!=nt:
msgs.append('tiled positions are not unique')
#end if
if len(msgs)>0:
msg = 'tiling check failed'
for m in msgs:
msg += '\n'+m
#end for
if exit:
self.error(msg)
#end if
#end if
return msg
#end def check_tiling
# test needed
[docs]
def kfold(self,tiling,kpoints,kweights):
if isinstance(tiling,int):
tiling = self.dim*[tiling]
#end if
tiling = np.array(tiling)
if tiling.shape==(self.dim,self.dim):
tiling = tiling.T
#end if
tilematrix,tilevector = reduce_tilematrix(tiling)
ncells = int(np.round( np.abs(det(tilematrix)) ))
kp = self.tile_points(kpoints,self.kaxes,tilematrix,tilevector)
kw = np.array(ncells*list(kweights),dtype=float)/ncells
return kp,kw
#end def kfold
[docs]
def get_smallest(self):
if self.has_folded():
return self.folded_structure
else:
return self
#end if
#end def get_smallest
# test needed
[docs]
def fold(self,small,*requests):
self.error('fold needs a developers attention to make it equivalent with tile')
if self.dim!=3:
self.error('fold is currently only implemented for 3 dimensions')
#end if
self.recenter_k()
corners = []
ndim = len(small.axes)
imin = np.empty((ndim,),dtype=int)
imax = np.empty((ndim,),dtype=int)
imin[:] = 1000000
imax[:] = -1000000
axinv = inv(self.kaxes)
center = self.kaxes.sum(0)/2
c = np.empty((1,3))
for k in -1,2:
for j in -1,2:
for i in -1,2:
c[:] = i,j,k
c = dot(c,small.kaxes)
u = dot(c-center,axinv)
for d in range(ndim):
imin[d] = min(int(np.floor(u[0,d])),imin[d])
imax[d] = max(int(np.ceil(u[0,d])),imax[d])
#end for
#end for
#end for
#end for
axes = small.kaxes
axinv = inv(small.kaxes)
center = small.kaxes.sum(0)/2
nkpoints = len(self.kpoints)
kindices = []
kpoints = []
shift = np.empty((ndim,))
kr = list(range(nkpoints))
for k in range(imin[2],imax[2]+1):
for j in range(imin[1],imax[1]+1):
for i in range(imin[0],imax[0]+1):
for n in kr:
shift[:] = i,j,k
shift = dot(shift,self.kaxes)
kp = self.kpoints[n]+shift
u = dot(kp-center,axinv)
if np.abs(u).max()<.5+1e-10:
kindices.append(n)
kpoints.append(kp)
#end if
#end for
#end for
#end for
#end for
kindices = np.array(kindices)
kpoints = np.array(kpoints)
inside = self.inside(kpoints,axes,center)
kindices = kindices[inside]
kpoints = kpoints[inside]
small.kpoints = kpoints
small.recenter_k()
kpoints = np.array(small.kpoints)
if len(requests)>0:
results = []
for request in requests:
if request=='kmap':
kmap = obj()
for k in self.kpoints:
kmap[tuple(k)] = []
#end for
for i in range(len(kpoints)):
kp = tuple(self.kpoints[kindices[i]])
kmap[kp].append(np.array(kpoints[i]))
#end for
for kl,ks in kmap.items():
kmap[kl] = np.array(ks)
#end for
res = kmap
elif request=='tilematrix':
res = self.tilematrix(small)
else:
self.error(request+' is not a recognized input to fold')
#end if
results.append(res)
#end if
return results
#end if
#end def fold
[docs]
def tilematrix(self,small=None,tol=1e-6,status=False):
if small is None:
if self.folded_structure is not None:
small = self.folded_structure
else:
return np.identity(self.dim,dtype=int)
#end if
#end if
tm = dot(self.axes,inv(small.axes))
tilemat = np.array(np.around(tm),dtype=int)
error = np.abs(tilemat-tm).sum()
non_integer_elements = error > tol
if status:
return tilemat,not non_integer_elements
else:
if non_integer_elements:
self.error(
"large cell cannot be constructed as an integer tiling of the small cell\n"
"large cell axes:\n"
+str(self.axes)+"\n"
"small cell axes: \n"
+str(small.axes)+"\n"
"large/small:\n"
+str(self.axes/small.axes)+"\n"
"tiling matrix:\n"
+str(tm)+"\n"
"integerized tiling matrix:\n"
+str(tilemat)+"\n"
"error: "+str(error)+"\n"
"tolerance: "+str(tol)
)
#end if
return tilemat
#end if
#end def tilematrix
[docs]
def primitive(self,source=None,tmatrix=False,add_kpath=False,**kwargs):
res = None
allowed_sources = set(['seekpath'])
if source is None or isinstance(source,bool):
source = 'seekpath'
#end if
if source not in allowed_sources:
self.error(
"source used to obtain primitive cell is unrecognized\n"
"source requested: {0}\n"
"allowed sources: {1}".format(
source, sorted(allowed_sources)
)
)
#end if
if source=='seekpath':
res_skp = get_seekpath_full(structure=self,primitive=True,**kwargs)
prim = res_skp.primitive
T = res_skp.prim_tmatrix
if add_kpath:
prim.add_kpoints(res_skp.explicit_kpoints_abs)
#end if
if tmatrix:
res = prim,T
else:
res = prim
#end if
else:
self.error(
'primitive source "{0}" is not implemented\n'
'please contact a developer'.format(source)
)
#end if
if prim.units!=self.units:
prim.change_units(self.units)
#end if
return res
#end def primitive
[docs]
def become_primitive(self,source=None,add_kpath=False,**kwargs):
prim = self.primitive(source=source,add_kpath=add_kpath,**kwargs)
self.clone_from(prim)
#end def become_primitive
[docs]
def add_kpoints(self,kpoints,kweights=None,unique=False,recenter=True,cell_unit=False):
if kweights is None:
kweights = np.ones((len(kpoints),))
#end if
if cell_unit:
kpoints = np.dot(np.array(kpoints),self.kaxes)
#end if
self.kpoints = np.append(self.kpoints,kpoints,axis=0)
self.kweights = np.append(self.kweights,kweights)
if unique:
self.unique_kpoints()
#end if
if recenter:
self.recenter_k() #added because qmcpack cannot handle kpoints outside the box
#end if
if self.is_tiled():
kp,kw = self.kfold(self.tmatrix,kpoints,kweights)
self.folded_structure.add_kpoints(kp,kw,unique=unique)
#end if
#end def add_kpoints
# test needed
[docs]
def clear_kpoints(self):
self.kpoints = np.empty((0,self.dim))
self.kweights = np.empty((0,))
if self.folded_structure is not None:
self.folded_structure.clear_kpoints()
#end if
#end def clear_kpoints
[docs]
def kgrid_from_kspacing(self,kspacing):
kgrid = []
for ka in self.kaxes:
km = norm(ka)
kg = int(np.ceil(km/kspacing))
kgrid.append(kg)
#end for
return tuple(kgrid)
#end def kgrid_from_kspacing
[docs]
def add_kmesh(self,kgrid=None,kshift=None,unique=False,kspacing=None):
if kspacing is not None:
kgrid = self.kgrid_from_kspacing(kspacing)
elif kgrid is None:
self.error('kgrid input is required by add_kmesh')
#end if
self.add_kpoints(kmesh(self.kaxes,kgrid,kshift),unique=unique)
#end def add_kmesh
[docs]
def add_symmetrized_kmesh(self,kgrid=None,kshift=(0,0,0),kspacing=None):
# find kgrid from kspacing, if requested
if kspacing is not None:
kgrid = self.kgrid_from_kspacing(kspacing)
elif kgrid is None:
self.error('kgrid input is required by add_kmesh')
#end if
# get spglib cell data structure
cell = self.spglib_cell()
# get the symmetry mapping
kmap,kpoints_int = spglib.get_ir_reciprocal_mesh(
kgrid,
cell,
is_shift=kshift
)
# create the Monkhorst-Pack mesh
kshift = np.array(kshift,dtype=float)
okgrid = 1.0/np.array(kgrid,dtype=float)
kpoints = np.empty(kpoints_int.shape,dtype=float)
for i,ki in enumerate(kpoints_int):
kpoints[i] = (ki+kshift)*okgrid
#end for
kpoints = dot(kpoints,self.kaxes)
# reduce to only the symmetric kpoints with weights
kwmap = obj()
for ik in kmap:
if ik not in kwmap:
kwmap[ik] = 1
else:
kwmap[ik] += 1
#end if
#end for
nkpoints = len(kwmap)
kpoints_symm = np.empty((nkpoints,self.dim),dtype=float)
kweights_symm = np.empty((nkpoints,),dtype=float)
n = 0
for ik,kw in kwmap.items():
kpoints_symm[n] = kpoints[ik]
kweights_symm[n] = kw
n+=1
#end for
self.add_kpoints(kpoints_symm,kweights_symm)
#end def add_symmetrized_kmesh
[docs]
def kpoints_unit(self,kpoints=None):
if kpoints is None:
kpoints = self.kpoints
#end if
return dot(kpoints,inv(self.kaxes))
#end def kpoints_unit
[docs]
def kpoints_reduced(self,kpoints=None):
if kpoints is None:
kpoints = self.kpoints
#end if
return kpoints*self.scale/(2*pi)
#end def kpoints_reduced
[docs]
def kpoints_qmcpack(self,kpoints=None):
if kpoints is None:
kpoints = self.kpoints.copy()
#end if
kpoints = self.recenter_k(kpoints,kcenter=(0,0,0))
kpoints = self.kpoints_unit(kpoints)
kpoints = -kpoints
return kpoints
#end def kpoints_qmcpack
# test needed
[docs]
def inversion_symmetrize_kpoints(self,tol=1e-10,folded=False):
kp = self.kpoints
kaxes = self.kaxes
ntable,dtable = self.neighbor_table(kp,-kp,kaxes,distances=True)
pairs = set()
keep = np.empty((len(kp),),dtype=bool)
keep[:] = True
for i in range(len(dtable)):
if keep[i] and dtable[i,0]<tol:
j = ntable[i,0]
if j!=i and keep[j]:
keep[j] = False
self.kweights[i] += self.kweights[j]
#end if
#end if
#end for
self.kpoints = self.kpoints[keep]
self.kweights = self.kweights[keep]
if folded and self.folded_structure is not None:
self.folded_structure.inversion_symmetrize_kpoints(tol)
#end if
#end def inversion_symmetrize_kpoints
# test needed
[docs]
def unique_points(self,points,axes,weights=None,tol=1e-10):
pmap = obj()
npoints = len(points)
if npoints>0:
if weights is None:
weights = np.ones((npoints,),dtype=int)
#end if
ntable,dtable = self.neighbor_table(points,points,axes,distances=True)
keep = np.empty((npoints,),dtype=bool)
keep[:] = True
pmo = obj()
for i in range(npoints):
if keep[i]:
pm = []
jn=0
while jn<npoints and dtable[i,jn]<tol:
j = ntable[i,jn]
pm.append(j)
if j!=i and keep[j]:
keep[j] = False
weights[i] += weights[j]
#end if
jn+=1
#end while
pmo[i] = set(pm)
#end if
#end for
points = points[keep]
weights = weights[keep]
j=0
for i in range(len(keep)):
if keep[i]:
pmap[j] = pmo[i]
j+=1
#end if
#end for
#end if
return points,weights,pmap
#end def unique_points
# test needed
[docs]
def unique_points_fast(self,points,axes,weights=None,tol=1e-10):
# use an O(N) cell table instead of an O(N^2) neighbor table
pmap = obj()
points = np.array(points)
axes = np.array(axes)
npoints = len(points)
if npoints>0:
if weights is None:
weights = np.ones((npoints,),dtype=int)
else:
weights = np.array(weights)
#end if
keep = np.ones((npoints,),dtype=bool)
# place all the points in the box, converted to unit coords
upoints = np.array(points)
axinv = inv(axes)
for i in range(len(points)):
u = dot(points[i],axinv)
upoints[i] = u-np.floor(u)
#end for
# create an integer array of cell indices
axmax = -1.0
for a in axes:
axmax = max(axmax,norm(a))
#end for
# make an integer space corresponding to 1e-7 self.units spatial resolution
cmax = np.uint64(1e7)*np.uint64(np.ceil(axmax))
ipoints = np.array(np.around(cmax*upoints),dtype=np.uint64)
ipoints[ipoints==cmax] = 0 # make the outer boundary the same as the inner boundary
# load the cell table with point indices
# points in the same cell are identical
ctable = obj()
i=0
for ip in ipoints:
ip = tuple(ip)
if ip not in ctable:
ctable[ip] = i
pmap[i] = [i]
else:
j = ctable[ip]
keep[i] = False
weights[j] += weights[i]
pmap[j].append(i)
#end if
i+=1
#end for
points = points[keep]
weights = weights[keep]
#end if
return points,weights,pmap
#end def unique_points_fast
# test needed
[docs]
def unique_positions(self,tol=1e-10,folded=False):
pos,weights,pmap = self.unique_points(self.pos,self.axes)
if len(pos)!=len(self.pos):
self.pos = pos
#end if
if folded and self.folded_structure is not None:
self.folded_structure.unique_positions(tol)
#end if
return pmap
#end def unique_positions
# test needed
[docs]
def unique_kpoints(self,tol=1e-10,folded=False):
kmap = obj()
kp = self.kpoints
if len(kp)>0:
kaxes = self.kaxes
ntable,dtable = self.neighbor_table(kp,kp,kaxes,distances=True)
npoints = len(kp)
keep = np.empty((len(kp),),dtype=bool)
keep[:] = True
kmo = obj()
for i in range(npoints):
if keep[i]:
km = []
jn=0
while jn<npoints and dtable[i,jn]<tol:
j = ntable[i,jn]
km.append(j)
if j!=i and keep[j]:
keep[j] = False
self.kweights[i] += self.kweights[j]
#end if
jn+=1
#end while
kmo[i] = set(km)
#end if
#end for
self.kpoints = self.kpoints[keep]
self.kweights = self.kweights[keep]
j=0
for i in range(len(keep)):
if keep[i]:
kmap[j] = kmo[i]
j+=1
#end if
#end for
#end if
if folded and self.folded_structure is not None:
self.folded_structure.unique_kpoints(tol)
#end if
return kmap
#end def unique_kpoints
[docs]
def kmap(self):
kmap = None
if self.folded_structure is not None:
fs = self.folded_structure
self.kpoints = np.array(fs.kpoints)
self.kweights = np.array(fs.kweights)
kmap = self.unique_kpoints()
#end if
return kmap
#end def kmap
# test needed
[docs]
def select_twist(self,selector='smallest',tol=1e-6):
index = None
invalid_selector = False
if isinstance(selector,str):
if selector=='smallest':
index = (self.kpoints**2).sum(1).argmin()
elif selector=='random':
index = randint(0,len(self.kpoints)-1)
else:
invalid_selector = True
#end if
elif isinstance(selector,(tuple,list,np.ndarray)):
ku_sel = np.array(selector,dtype=float)
n = 0
for ku in self.kpoints_unit():
if norm(ku-ku_sel)<tol:
index = n
break
#end if
n+=1
#end for
if index is None:
self.error(
"cannot identify twist number\n"
"twist requested: {0}\n"
"twists present: {1}".format(
ku_sel, sorted([tuple(k) for k in self.kpoints_unit()])
)
)
#end if
else:
invalid_selector = True
#end if
if invalid_selector:
self.error(
"cannot identify twist number\n"
"invalid selector provided: {0}\n"
"valid string inputs for selector: smallest, random\n"
"selector can also be a length 3 tuple, list or array (a twist vector)".format(selector)
)
#end if
return index
#end def select_twist
# test needed
[docs]
def fold_pos(self,large,tol=0.001):
vratio = large.volume()/self.volume()
if np.abs(vratio-int(np.around(vratio)))>1e-6:
self.error(
"cannot fold positions from large cell into current one\n"
"large cell volume is not an integer multiple of the current one\n"
"large cell volume: {0}\n"
"current cell volume: {1}\n"
"volume ratio: {2}".format(large.volume(),self.volume(),vratio)
)
T,success = large.tilematrix(self,status=True)
if not success:
self.error(
"cannot fold positions from large cell into current one\n"
"cells are related by non-integer tilematrix"
)
#end if
nnearest = int(np.around(vratio))
self.elem = large.elem.copy()
self.pos = large.pos.copy()
self.recenter()
nt,dt = self.neighbor_table(distances=True)
nt = nt[:,:nnearest]
dt = dt[:,:nnearest]
if dt.ravel().max()>tol:
self.error(
"cannot fold positions from large cell into current one\n"
"positions of equivalent atoms are further apart than the tolerance\n"
"max distance encountered: {0}\n"
"tolerance: {1}".format(dt.ravel().max(),tol)
)
#end if
counts = np.zeros((len(self.pos),),dtype=int)
for n in nt.ravel():
counts[n] += 1
#end for
if (counts!=nnearest).any():
self.error(
"cannot fold positions from large cell into current one\n"
"each atom must have {0} equivalent positions\n"
"some atoms found with the following equivalent position counts: {1}".format(
nnearest, counts[counts!=nnearest]
)
)
#end if
ind_visited = set()
neigh_map = obj()
keep = []
n=0
for nset in nt:
if n not in ind_visited:
neigh_map[n] = nset
keep.append(n)
for ind in nset:
ind_visited.add(ind)
#end for
#end if
n+=1
#end for
if len(ind_visited)!=len(self.pos):
self.error(
"cannot fold positions from large cell into current one\n"
"some equivalent atoms could not be identified"
)
#end if
new_elem = []
new_pos = []
for n in keep:
nset = neigh_map[n]
elist = list(set(self.elem[nset]))
if len(elist)!=1:
self.error(
"cannot fold positions from large cell into current one\n"
"species of some equivalent atoms do not match"
)
#end if
new_elem.append(elist[0])
new_pos.append(self.pos[nset].mean(0))
#end for
self.set_elem(new_elem)
self.set_pos(new_pos)
#end def fold_pos
[docs]
def pos_unit(self,pos=None):
if pos is None:
pos = self.pos
#end if
return dot(pos,inv(self.axes))
#end def pos_unit
[docs]
def pos_to_cartesian(self):
self.pos = dot(self.pos,self.axes)
#end def pos_to_cartesian
# test needed
[docs]
def at_Gpoint(self):
kpu = self.kpoints_unit()
kg = np.array([0,0,0])
return len(kpu)==1 and norm(kg-kpu[0])<1e-6
#end def at_Gpoint
# test needed
[docs]
def at_Lpoint(self):
kpu = self.kpoints_unit()
kg = np.array([.5,.5,.5])
return len(kpu)==1 and norm(kg-kpu[0])<1e-6
#end def at_Lpoint
# test needed
[docs]
def at_real_kpoint(self):
kpu = 2*self.kpoints_unit()
return len(kpu)==1 and np.abs(kpu-np.around(kpu)).sum()<1e-6
#end def at_real_kpoint
# test needed
[docs]
def bonds(self,neighbors,vectors=False):
if self.dim!=3:
self.error('bonds is currently only implemented for 3 dimensions')
#end if
natoms,dim = self.pos.shape
centers = np.empty((natoms,neighbors,dim))
distances = np.empty((natoms,neighbors))
vect = np.empty((natoms,neighbors,dim))
t = self.tile((3,3,3))
t.recenter(self.center)
nn = nearest_neighbors(neighbors+1,t.pos,self.pos)
for i in range(natoms):
ii = nn[i,0]
n=0
for jj in nn[i,1:]:
p1 = t.pos[ii]
p2 = t.pos[jj]
centers[i,n,:] = (p1+p2)/2
distances[i,n]= sqrt(((p1-p2)**2).sum())
vect[i,n,:] = p2-p1
n+=1
#end for
#end for
sn = self.copy()
nnr = nn[:,1:].ravel()
sn.set_elem(t.elem[nnr])
sn.pos = t.pos[nnr]
sn.recenter()
indices = self.locate(sn.pos)
indices = indices.reshape(natoms,neighbors)
if not vectors:
return indices,centers,distances
else:
return indices,centers,distances,vect
#end if
#end def bonds
# test needed
[docs]
def displacement(self,reference,map=False):
if self.dim!=3:
self.error('displacement is currently only implemented for 3 dimensions')
#end if
ref = reference.tile((3,3,3))
ref.recenter(reference.center)
rmap = np.array(3**3*list(range(len(reference.pos)),dtype=int))
nn = nearest_neighbors(1,ref.pos,self.pos).ravel()
displacement = self.pos - ref.pos[nn]
if not map:
return displacement
else:
return displacement,rmap[nn]
#end if
#end def displacement
# test needed
[docs]
def scalar_displacement(self,reference):
return sqrt((self.displacement(reference)**2).sum(1))
#end def scalar_displacement
# test needed
[docs]
def distortion(self,reference,neighbors):
if self.dim!=3:
self.error('distortion is currently only implemented for 3 dimensions')
#end if
if reference.volume()/self.volume() < 1.1:
ref = reference.tile((3,3,3))
ref.recenter(reference.center)
else:
ref = reference
#end if
rbi,rbc,rbl,rbv = ref.bonds(neighbors,vectors=True)
sbi,sbc,sbl,sbv = self.bonds(neighbors,vectors=True)
nn = nearest_neighbors(1,reference.pos,self.pos).ravel()
distortion = np.empty(sbv.shape)
magnitude = np.empty((len(self.pos),))
for i in range(len(self.pos)):
ir = nn[i]
bonds = sbv[i]
rbonds = rbv[ir]
ib = np.empty((neighbors,),dtype=int)
ibr = np.empty((neighbors,),dtype=int)
r = list(range(neighbors))
rr = list(range(neighbors))
for n in range(neighbors):
mindist = 1e99
ibmin = -1
ibrmin = -1
for nb in r:
for nbr in rr:
d = norm(bonds[nb]-rbonds[nbr])
if d<mindist:
mindist=d
ibmin=nb
ibrmin=nbr
#end if
#end for
#end for
ib[n]=ibmin
ibr[n]=ibrmin
r.remove(ibmin)
rr.remove(ibrmin)
#end for
#end for
d = bonds[ib]-rbonds[ibr]
distortion[i] = d
magnitude[i] = (sqrt((d**2).sum(axis=1))).sum()
#end for
return distortion,magnitude
#end def distortion
# test needed
[docs]
def bond_compression(self,reference,neighbors):
ref = reference
rbi,rbc,rbl = ref.bonds(neighbors)
sbi,sbc,sbl = self.bonds(neighbors)
bondlen = rbl.mean()
return np.abs(1.-sbl/bondlen).max(axis=1)
#end def bond_compression
# test needed
[docs]
def boundary(self,dims=(0,1,2),dtol=1e-6):
dim_eff = len(dims)
natoms,dim = self.pos.shape
bdims = np.array(dim*[False])
for d in dims:
bdims[d] = True
#end for
p = self.pos[:,bdims]
indices = convex_hull(p,dim_eff,dtol)
return indices
#end def boundary
[docs]
def embed(self,small,dims=(0,1,2),dtol=1e-6,utol=1e-6):
small = small.copy()
small.recenter()
center = np.array(self.center)
self.recenter(small.center)
bind = small.boundary(dims,dtol)
bpos = small.pos[bind]
belem= small.elem[bind]
nn = nearest_neighbors(1,self.pos,bpos).ravel()
mpos = self.pos[nn]
dr = (mpos-bpos).mean(0)
for i in range(len(bpos)):
bpos[i]+=dr
#end for
dmax = sqrt(((mpos-bpos)**2).sum(1)).max()
for i in range(len(small.pos)):
small.pos[i]+=dr
#end for
ins,surface = small.inside(self.pos,tol=utol,separate=True)
replaced = np.empty((len(self.pos),),dtype=bool)
replaced[:] = False
inside = replaced.copy()
inside[ins] = True
nn = nearest_neighbors(1,self.pos,small.pos).ravel()
elist = list(self.elem)
plist = list(self.pos)
pos = small.pos
elem = small.elem
for i in range(len(pos)):
n = nn[i]
if not replaced[n]:
elist[n] = elem[i]
plist[n] = pos[i]
replaced[n] = True
else:
elist.append(elem[i])
plist.append(pos[i])
#end if
#end for
remove = np.arange(len(self.pos))[inside & np.logical_not(replaced)]
remove.sort()
remove = np.flipud(remove)
for i in remove:
elist.pop(i)
plist.pop(i)
#end for
self.set_elem(elist)
self.pos = np.array(plist)
self.recenter(center)
return dmax
#end def embed
# test needed
[docs]
def shell(self,cell,neighbors,direction='in'):
if self.dim!=3:
self.error('shell is currently only implemented for 3 dimensions')
#end if
dd = {'in':equate,'out':negate}
dir = dd[direction]
natoms,dim=self.pos.shape
ncells=3**3
ntile = ncells*natoms
pos = np.empty((ntile,dim))
ind = np.empty((ntile,),dtype=int)
oind = list(range(natoms))
for nt in range(ncells):
n=nt*natoms
ind[n:n+natoms]=oind[:]
pos[n:n+natoms]=self.pos[:]
#end for
nt=0
for k in -1,0,1:
for j in -1,0,1:
for i in -1,0,1:
iv = np.array([[i,j,k]])
v = dot(iv,self.axes)
for d in range(dim):
ns = nt*natoms
ne = ns+natoms
pos[ns:ne,d] += v[0,d]
#end for
nt+=1
#end for
#end for
#end for
inside = np.empty(ntile,)
inside[:]=False
ins = cell.inside(pos)
inside[ins]=True
iishell = set()
nn = nearest_neighbors(neighbors,pos)
for ii in range(len(nn)):
for jj in nn[ii]:
in1 = inside[ii]
in2 = inside[jj]
if dir(in1 and not in2):
iishell.add(ii)
#end if
if dir(in2 and not in1):
iishell.add(jj)
#end if
#end if
#end if
ishell = ind[list(iishell)]
return ishell
#end def shell
[docs]
def interpolate(self,other,images,min_image=True,recenter=True,match_com=False,chained=False):
s1 = self.copy()
s2 = other.copy()
s1.remove_folded()
s2.remove_folded()
if s2.units!=s1.units:
s2.change_units(s1.units)
#end if
if (s1.elem!=s2.elem).any():
self.error(
"cannot interpolate structures, atoms do not match\n"
" atoms1: {0}\n"
" atoms2: {1}".format(
s1.elem, s2.elem
)
)
#end if
structures = []
npath = images+2
c1 = s1.center
c2 = s2.center
ax1 = s1.axes
ax2 = s2.axes
pos1 = s1.pos
pos2 = s2.pos
min_image &= np.abs(ax1-ax2).max()<1e-6
if min_image:
dp = self.min_image_vectors(pos1,pos2,ax1,pairs=False)
pos2 = pos1 + dp
#end if
if match_com:
com1 = pos1.mean(axis=0)
com2 = pos2.mean(axis=1)
dcom = com1-com2
for n in range(len(pos2)):
pos2[n] += dcom
#end for
if chained:
other.pos = pos2
#end if
#end if
for n in range(npath):
f1 = 1.-float(n)/(npath-1)
f2 = 1.-f1
center = f1*c1 + f2*c2
axes = f1*ax1 + f2*ax2
pos = f1*pos1 + f2*pos2
s = s1.copy()
s.reset_axes(axes)
s.center = center
s.pos = pos
if recenter:
s.recenter()
#end if
structures.append(s)
#end for
return structures
#end def interpolate
# returns madelung potential constant v_M
# see equation 7 in PRB 78 125106 (2008)
[docs]
def madelung(self,axes=None,tol=1e-10):
if self.dim!=3:
self.error('madelung is currently only implemented for 3 dimensions')
#end if
if axes is None:
a = self.axes.T.copy()
else:
a = axes.T.copy()
#end if
if self.units!='B':
a = convert(a,self.units,'B')
#end if
volume = np.abs(det(a))
b = 2*pi*inv(a).T
rconv = 8*(3.*volume/(4*pi))**(1./3)
kconv = 2*pi/rconv
gconst = -1./(4*kconv**2)
vmc = -pi/(kconv**2*volume)-2*kconv/sqrt(pi)
nshells = 20
vshell = [0.]
p = Sobj()
m = Sobj()
for n in range(1,nshells+1):
i = np.mgrid[-n:n+1,-n:n+1,-n:n+1]
i = i.reshape(3,(2*n+1)**3)
R = sqrt((dot(a,i)**2).sum(0))
G2 = (dot(b,i)**2).sum(0)
izero = n + n*(2*n+1) + n*(2*n+1)**2
p.R = R[0:izero]
p.G2 = G2[0:izero]
m.R = R[izero+1:]
m.G2 = G2[izero+1:]
domains = [p,m]
vshell.append(0.)
for d in domains:
vshell[n] += (erfc(kconv*d.R)/d.R).sum() + 4*pi/volume*(np.exp(gconst*d.G2)/d.G2).sum()
#end for
if np.abs(vshell[n]-vshell[n-1])<tol:
break
#end if
#end for
vm = vmc + vshell[-1]
if axes is None:
self.Vmadelung = vm
#end if
return vm
#end def madelung
[docs]
def makov_payne(self,q=1,eps=1.0,units='Ha',order=1):
if order!=1:
self.error('Only first order Makov-Payne correction is currently supported.')
#end if
if 'Vmadelung' not in self:
vm = self.madelung()
else:
vm = self.Vmadelung
#end if
mp = -0.5*q**2*vm/eps
if units!='Ha':
mp = convert(mp,'Ha',units)
#end if
return mp
#end def makov_payne
[docs]
def read(self,filepath,format=None,elem=None,block=None,grammar='1.1',cell='prim',contents=False):
if os.path.exists(filepath):
filepath = path_string(filepath)
path,file = os.path.split(filepath)
if format is None:
if '.' in file:
name,format = file.rsplit('.',1)
elif file.lower().endswith('poscar'):
format = 'poscar'
else:
self.error(
"file format could not be determined\n"
"unrecognized file: {0}".format(filepath)
)
#end if
#end if
elif not contents:
self.error('file does not exist: {0}'.format(filepath))
#end if
if format is None:
self.error('file format must be provided')
#end if
self.mag = None
self.frozen = None
format = format.lower()
if format=='xyz':
self.read_xyz(filepath)
elif format=='xsf':
self.read_xsf(filepath)
elif format=='poscar':
self.read_poscar(filepath,elem=elem)
elif format=='cif':
self.read_cif(filepath,block=block,grammar=grammar,cell=cell)
elif format=='fhi-aims':
self.read_fhi_aims(filepath)
else:
self.error(
"cannot read structure from file\n"
"unsupported file format: {0}".format(format)
)
#end if
if self.has_axes():
self.set_bconds('ppp')
#end if
#end def read
[docs]
def read_xyz(self,filepath):
elem = []
pos = []
if os.path.exists(filepath):
with open(filepath, "r") as f:
lines = f.read().strip().splitlines()
else:
lines = filepath.strip().splitlines() # "filepath" is file contents
#end if
if len(lines)>1:
ntot = int(lines[0].strip())
natoms = 0
e = None
p = None
try:
tokens = lines[1].split()
if len(tokens)==4:
e = tokens[0]
p = np.array(tokens[1:],float)
#end if
except:
pass
#end try
if p is not None:
elem.append(e)
pos.append(p)
natoms+=1
#end if
if len(lines)>2:
for l in lines[2:]:
tokens = l.split()
if len(tokens)==4:
elem.append(tokens[0])
pos.append(np.array(tokens[1:],float))
natoms+=1
if natoms==ntot:
break
#end if
#end if
#end for
#end if
if natoms!=ntot:
self.error(
"xyz file read failed\n"
"attempted to read file: {0}\n"
"number of atoms expected: {1}\n"
"number of atoms found: {2}".format(
filepath, ntot, natoms
)
)
#end if
#end if
self.dim = 3
self.set_elem(elem)
self.pos = np.array(pos)
self.units = 'A'
#end def read_xyz
[docs]
def read_xsf(self,filepath):
if isinstance(filepath,XsfFile):
f = filepath
elif os.path.exists(filepath):
f = XsfFile(filepath)
else:
f = XsfFile()
f.read_text(filepath) # "filepath" is file contents
#end if
elem = []
for n in f.elem:
if isinstance(n,str):
elem.append(n)
else:
elem.append(Elements(n).symbol)
#end if
#end for
self.dim = 3
self.units = 'A'
self.reset_axes(f.primvec)
self.set_elem(elem)
self.pos = f.pos
#end def read_xsf
[docs]
def read_poscar(self,filepath,elem=None):
if os.path.exists(filepath):
with open(filepath, "r") as f:
lines = f.read().splitlines()
else:
lines = filepath.splitlines() # "filepath" is file contents
#end if
nlines = len(lines)
min_lines = 8
if nlines<min_lines:
self.error(
"POSCAR file must have at least {0} lines\n"
" only {1} lines found".format(min_lines, nlines)
)
#end if
dim = 3
scale = float(lines[1].strip())
axes = np.empty((dim,dim))
axes[0] = np.array(lines[2].split(),dtype=float)
axes[1] = np.array(lines[3].split(),dtype=float)
axes[2] = np.array(lines[4].split(),dtype=float)
if scale<0.0:
scale = np.abs(scale)/det(axes)
#end if
axes = scale*axes
tokens = lines[5].split()
if tokens[0].isdigit():
counts = np.array(tokens,dtype=int)
if elem is None:
self.error(
"variable elem must be provided to read_poscar() to assign atomic species to positions for POSCAR format"
)
elif len(elem)!=len(counts):
self.error(
"one elem must be given for each element count in the POSCAR file\n"
" number of elem counts: {0}\n"
" number of elem given: {1}".format(
len(counts), len(elem)
)
)
#end if
lcur = 6
else:
elem = tokens
counts = np.array(lines[6].split(),dtype=int)
lcur = 7
#end if
species = elem
# relabel species that have multiple occurances
sset = set(species)
for spec in sset:
if species.count(spec)>1:
cnt=0
for n in range(len(species)):
specn = species[n]
if specn==spec:
cnt+=1
species[n] = specn+str(cnt)
#end if
#end for
#end if
#end for
elem = []
for i in range(len(counts)):
elem.extend(counts[i]*[species[i]])
#end for
self.dim = dim
self.units = 'A'
self.reset_axes(axes)
if lcur<len(lines) and len(lines[lcur])>0:
c = lines[lcur].lower().strip()[0]
lcur+=1
else:
return
#end if
selective_dynamics = c=='s'
if selective_dynamics: # Selective dynamics
if lcur<len(lines) and len(lines[lcur])>0:
c = lines[lcur].lower().strip()[0]
lcur+=1
else:
return
#end if
#end if
cartesian = c=='c' or c=='k'
npos = counts.sum()
if lcur+npos>len(lines):
return
#end if
spos = []
for i in range(npos):
spos.append(lines[lcur+i].split())
#end for
spos = np.array(spos)
pos = np.array(spos[:,0:3],dtype=float)
if cartesian:
pos = scale*pos
else:
pos = dot(pos,axes)
#end if
self.set_elem(elem)
self.pos = pos
if selective_dynamics or spos.shape[1]>3:
move = np.array(spos[:,3:6],dtype=str)
self.freeze(list(range(self.size())),directions=move=='F')
#end if
#end def read_poscar
[docs]
def read_cif(self,filepath,block=None,grammar='1.1',cell='prim'):
axes,elem,pos,units = read_cif(filepath,block,grammar,cell,args_only=True)
self.dim = 3
self.set_axes(axes)
self.set_elem(elem)
self.pos = pos
self.units = units
#end def read_cif
# test needed
[docs]
def read_fhi_aims(self,filepath):
if os.path.exists(filepath):
with open(filepath, "r") as f:
lines = f.read().splitlines()
else:
lines = filepath.splitlines() # "filepath" is contents
#end if
axes = []
pos = []
elem = []
constrain_relax = []
unit_pos = False
for line in lines:
ls = line.strip()
if len(ls)>0 and ls[0]!='#':
tokens = ls.split()
t0 = tokens[0]
if t0=='lattice_vector':
axes.append(tokens[1:])
elif t0=='atom_frac':
pos.append(tokens[1:4])
elem.append(tokens[4])
unit_pos = True
elif t0=='atom':
pos.append(tokens[1:4])
elem.append(tokens[4])
elif t0=='constrain_relaxation':
constrain_relax.append(tokens[1])
elif t0.startswith('initial'):
None
else:
#None
self.error(
"unrecogonized or not yet supported token in fhi-aims geometry file: {0}".format(t0)
)
#end if
#end if
#end for
axes = np.array(axes,dtype=float)
pos = np.array(pos,dtype=float)
if unit_pos:
pos = dot(pos,axes)
#end if
self.dim = 3
if len(axes)>0:
self.set_axes(axes)
#end if
self.set_elem(elem)
self.pos = pos
self.units = 'A'
if len(constrain_relax)>0:
constrain_relax = np.array(constrain_relax)
self.freeze(list(range(self.size())),directions=constrain_relax=='.true.')
#end if
#end def read_fhi_aims
[docs]
def write(self,filepath=None,format=None):
if filepath is None and format is None:
self.error('please specify either the filepath or format arguments to write()')
elif format is None:
if isinstance(filepath, Path):
format = filepath.suffix.lstrip(".")
elif '.' in filepath:
format = filepath.split('.')[-1]
else:
self.error(
"file format could not be determined\n"
"either request the format directly with the format keyword or add a file format extension to the file name"
)
#end if
#end if
format = format.lower()
if format=='xyz':
c = self.write_xyz(filepath)
elif format=='xsf':
c = self.write_xsf(filepath)
elif format=='poscar':
c = self.write_poscar(filepath)
elif format=='fhi-aims':
c = self.write_fhi_aims(filepath)
else:
self.error('file format {0} is unrecognized'.format(format))
#end if
return c
#end def write
[docs]
def pos_to_str(self, units: str = "A", with_elem: bool = False):
"""Write the positions of a structure to a string, optionally with atomic symbols.
Parameters
----------
units : str, default="A"
Units for the positions.
with_elem : bool, default=False
Optionally write the atomic symbols with the positions.
Notes
-----
This function will write the positions in the format
(in this case using ``with_elem=True``)
.. code-block:: python
"{element:2} {dim1:12.8f} {dim2:12.8f} ... {dimN:12.8f}\\n"
"""
s = self.copy()
s.change_units(units)
c = ""
for i in range(len(s.elem)):
e = s.elem[i]
p = s.pos[i]
if with_elem:
c += f"{e:2}"
for j in p:
c += f"{j:12.8f} "
c = c.rstrip()
c += "\n"
#end for
return c
#end def pos_to_str
[docs]
def write_xyz(self, filepath=None):
"""Write a ``Structure`` object to an XYZ file
Parameters
----------
filepath : PathLike or None, default=None
Path to where the XYZ file should be written.
If this is ``None``, then this function just returns
what would have been written to the XYZ file.
Returns
-------
xyz : str
The text that was or would have been written to the XYZ file.
Notes
-----
To get a string of only the atomic positions, use ``pos_to_str()`` instead.
"""
if self.dim!=3:
self.error('write_xyz is currently only implemented for 3 dimensions')
#end if
s = self.copy()
s.change_units("A")
c = ''
c += str(len(s.elem))+'\n\n'
for i in range(len(s.elem)):
e = s.elem[i]
p = s.pos[i]
c+=' {0:2} {1:12.8f} {2:12.8f} {3:12.8f}\n'.format(e,p[0],p[1],p[2])
#end for
if filepath is not None:
with open(filepath, "w") as f:
f.write(c)
#end if
return c
#end def write_xyz
[docs]
def write_xsf(self,filepath=None):
if self.dim!=3:
self.error('write_xsf is currently only implemented for 3 dimensions')
#end if
s = self.copy()
s.change_units('A')
c = ' CRYSTAL\n'
c += ' PRIMVEC\n'
for a in s.axes:
c += ' {0:12.8f} {1:12.8f} {2:12.8f}\n'.format(*a)
#end for
c += ' PRIMCOORD\n'
c += ' {0} 1\n'.format(len(s.elem))
for i in range(len(s.elem)):
e = s.elem[i]
identified, element = Elements.is_element(e, return_element=True)
if not identified:
self.error(
"{0} is not an element\n"
"xsf file cannot be written".format(e)
)
#end if
enum = element.atomic_number
r = s.pos[i]
c += ' {0:>3} {1:12.8f} {2:12.8f} {3:12.8f}\n'.format(enum,r[0],r[1],r[2])
#end for
if filepath is not None:
with open(filepath, "w") as f:
f.write(c)
#end if
return c
#end def write_xsf
[docs]
def write_poscar(self,filepath=None):
s = self.copy()
s.change_units('A')
species,species_count = s.order_by_species()
poscar = PoscarFile()
poscar.scale = 1.0
poscar.axes = s.axes
poscar.elem = species
poscar.elem_count = species_count
poscar.coord = 'cartesian'
poscar.pos = s.pos
c = poscar.write_text()
if filepath is not None:
with open(filepath, "w") as f:
f.write(c)
#end if
return c
#end def write_poscar
# test needed
[docs]
def write_fhi_aims(self,filepath=None):
s = self.copy()
s.change_units('A')
c = ''
c+='\n'
for a in s.axes:
c += 'lattice_vector {0: 12.8f} {1: 12.8f} {2: 12.8f}\n'.format(*a)
#end for
c+='\n'
for p,e in zip(self.pos,self.elem):
c += 'atom_frac {0: 12.8f} {1: 12.8f} {2: 12.8f} {3}\n'.format(p[0],p[1],p[2],e)
#end for
if filepath is not None:
with open(filepath, "w") as f:
f.write(c)
#end if
return c
#end def write_fhi_aims
[docs]
def plot2d_ax(self,ix=0,iy=1,*args,**kwargs):
if self.dim!=3:
self.error('plot2d_ax is currently only implemented for 3 dimensions')
#end if
iz = list(set([0,1,2])-set([ix,iy]))[0]
ax = self.axes.copy()
a = self.axes[iz]
dc = self.center-ax.sum(0)/2
pp = np.array([0*a,ax[ix],ax[ix]+ax[iy],ax[iy],0*a])
for i in range(len(pp)):
pp[i]+=dc
pp[i]-=dot(a,pp[i])/dot(a,a)*a
#end for
plot(pp[:,ix],pp[:,iy],*args,**kwargs)
#end def plot2d_ax
[docs]
def plot2d_pos(self,ix=0,iy=1,*args,**kwargs):
if self.dim!=3:
self.error('plot2d_pos is currently only implemented for 3 dimensions')
#end if
iz = list(set([0,1,2])-set([ix,iy]))[0]
pp = self.pos.copy()
a = self.axes[iz]
for i in range(len(pp)):
pp[i] -= dot(a,pp[i])/dot(a,a)*a
#end for
plot(pp[:,ix],pp[:,iy],*args,**kwargs)
#end def plot2d_pos
[docs]
def plot2d_points(self,points,ix=0,iy=1,*args,**kwargs):
if self.dim!=3:
self.error('plot2d_points is currently only implemented for 3 dimensions')
#end if
iz = list(set([0,1,2])-set([ix,iy]))[0]
pp = np.array(points,dtype=float)
a = self.axes[iz]
for i in range(len(pp)):
pp[i] -= dot(a,pp[i])/dot(a,a)*a
#end for
plot(pp[:,ix],pp[:,iy],*args,**kwargs)
#end def plot2d_points
[docs]
def plot2d(self,pos_style='b.',ax_style='k-'):
if self.dim!=3:
self.error('plot2d is currently only implemented for 3 dimensions')
#end if
subplot(1,3,1)
self.plot2d_ax(0,1,ax_style,lw=2)
self.plot2d_pos(0,1,pos_style)
title('a1,a2')
subplot(1,3,2)
self.plot2d_ax(1,2,ax_style,lw=2)
self.plot2d_pos(1,2,pos_style)
title('a2,a3')
subplot(1,3,3)
self.plot2d_ax(2,0,ax_style,lw=2)
self.plot2d_pos(2,0,pos_style)
title('a3,a1')
#end def plot2d
[docs]
def plot2d_kax(self,ix,iy,*args,**kwargs):
if self.dim!=3:
self.error('plot2d_ax is currently only implemented for 3 dimensions')
#end if
iz = list(set([0,1,2])-set([ix,iy]))[0]
ax = self.kaxes.copy()
a = ax[iz]
dc = 0*a
pp = np.array([0*a,ax[ix],ax[ix]+ax[iy],ax[iy],0*a])
for i in range(len(pp)):
pp[i]+=dc
pp[i]-=dot(a,pp[i])/dot(a,a)*a
#end for
plot(pp[:,ix],pp[:,iy],*args,**kwargs)
#end def plot2d_kax
[docs]
def plot2d_kp(self,ix,iy,*args,**kwargs):
if self.dim!=3:
self.error('plot2d_kp is currently only implemented for 3 dimensions')
#end if
iz = list(set([0,1,2])-set([ix,iy]))[0]
pp = self.kpoints.copy()
a = self.kaxes[iz]
for i in range(len(pp)):
pp[i] -= dot(a,pp[i])/dot(a,a)*a
#end for
plot(pp[:,ix],pp[:,iy],*args,**kwargs)
#end def plot2d_kp
[docs]
def show(self,viewer='vmd',filepath='/tmp/tmp.xyz'):
if self.dim!=3:
self.error('show is currently only implemented for 3 dimensions')
#end if
self.write_xyz(filepath)
os.system(viewer+' '+filepath)
#end def show
# minimal ASE Atoms-like interface to Structure objects for spglib
[docs]
def get_cell(self):
return self.axes.copy()
#end def get_cell
[docs]
def get_scaled_positions(self):
return self.pos_unit()
#end def get_scaled_positions
[docs]
def get_number_of_atoms(self):
return len(self.elem)
#end def get_number_of_atoms
[docs]
def get_atomic_numbers(self):
an = []
for e in self.elem:
iselem, element = Elements.is_element(e, return_element=True)
if not iselem:
self.error('Atomic symbol, {}, not recognized'.format(element))
else:
an.append(element.atomic_number)
#end if
#end for
return np.array(an,dtype='intc')
#end def get_atomic_numbers
[docs]
def get_magnetic_moments(self):
self.error('structure objects do not currently support magnetic moments')
#end def get_magnetic_moments
# direct spglib interface
[docs]
def spglib_cell(self):
lattice = self.axes.copy()
positions = self.pos_unit()
numbers = self.get_atomic_numbers()
cell = (lattice,positions,numbers)
return cell
#end def spglib_cell
[docs]
def get_symmetry(self,symprec=1e-5):
cell = self.spglib_cell()
return spglib.get_symmetry(cell,symprec=symprec)
#end def get_symmetry
[docs]
def get_symmetry_dataset(self,symprec=1e-5, angle_tolerance=-1.0, hall_number=0):
cell = self.spglib_cell()
ds = spglib.get_symmetry_dataset(
cell,
symprec = symprec,
angle_tolerance = angle_tolerance,
hall_number = hall_number,
)
return ds
#end def get_symmetry
# functions based on direct spglib interface
[docs]
def symmetry_data(self,*args,**kwargs):
ds = self.get_symmetry_dataset(*args,**kwargs)
if isinstance(ds,dict):
# Spglib version < v.2.5.0, see https://spglib.readthedocs.io/en/stable/releases.html
ds = obj(ds)
elif isinstance(ds,spglib.SpglibDataset):
# Spglib version >= v.2.5.0
ds = obj(ds.__dict__)
else:
raise TypeError('Invalid symmetry dataset type: {}'.format(type(ds)))
for k,v in ds.items():
if isinstance(v,dict):
ds[k] = obj(v)
#end if
#end for
return ds
#end def symmetry_data
[docs]
def bravais_lattice_name(self,symm_data=None):
if symm_data is None:
symm_data = self.symmetry_data()
#end if
sg = symm_data.number
name = symm_data.international
if not isinstance(sg,int) or sg<1 or sg>230:
self.error('Invalid space group from spglib: {}'.format(sg))
#end if
if not isinstance(name,str):
self.error('Invalid space group name from spglib: {}'.format(name))
#end if
bv = None
if sg>=1 and sg<=2:
bv = 'triclinic_'+name[0]
elif sg>=3 and sg<=15:
bv = 'monoclinic_'+name[0]
elif sg>=16 and sg<=74:
bv = 'orthorhombic_'+name[0]
elif sg>=75 and sg<=142:
bv = 'tetragonal_'+name[0]
elif sg>=143 and sg<=167:
bv = 'trigonal_'+name[0]
elif sg>=168 and sg<=194:
bv = 'hexagonal_'+name[0]
elif sg>=195 and sg<=230:
bv = 'cubic_'+name[0]
#end if
if bv is None:
self.error(
"Bravais lattice could not be determined.\n"
"Space group number and name: {} {}".format(sg, name)
)
#end if
return bv
#end def bravais_lattice_name
# test needed
[docs]
def space_group_operations(self,tol=1e-5,unit=False):
ds = self.get_symmetry(symprec=tol)
if ds is None:
self.error(
"Symmetry search failed.\n"
"spglib error message:\n"
"{}".format(spglib.get_error_message())
)
#end if
ds = obj(ds)
rotations = ds.rotations
translations = ds.translations
if not unit:
# Transform to Cartesian
axes = self.axes
axinv = inv(axes)
for n,(R,t) in enumerate(zip(rotations,translations)):
rotations[n] = np.dot(axinv,np.dot(R,axes))
translations[n] = np.dot(t,axes)
#end for
#end if
return rotations,translations
#end def space_group_operations
[docs]
def point_group_operations(self,tol=1e-5,unit=False):
rotations,translations = self.space_group_operations(tol=tol,unit=unit)
no_trans = translations.max(axis=1) < tol
return rotations[no_trans]
#end def point_group_operations
[docs]
def check_point_group_operations(self,rotations=None,tol=1e-5,unit=False,dtol=1e-5,ncheck=1,exit=False):
if rotations is None:
rotations = self.point_group_operations(tol=tol,unit=unit)
#ned if
r = self.pos
if ncheck=='all':
ncheck = len(r)
#end if
all_same = True
for n in range(ncheck):
rc = r[n]
for R in rotations:
rp = np.dot(r-rc,R)+rc
dt = self.min_image_distances(r,rp)
same = True
for d in dt:
same &= dt.min()<dtol
#end for
all_same &= same
#end for
#end for
if not all_same and exit:
self.error('Point group operators are not all symmetries of the structure.')
#end if
return all_same
#end def check_point_group_operations
[docs]
def equivalent_atoms(self):
ds = self.symmetry_data()
# collect sets of species labels
species_by_specnum = obj()
for e,sn in zip(self.elem,ds.equivalent_atoms):
is_elem,element = Elements.is_element(e, return_element=True)
if sn not in species_by_specnum:
species_by_specnum[sn] = set()
#end if
species_by_specnum[sn].add(element.symbol)
#end for
for sn,sset in species_by_specnum.items():
if len(sset)>1:
self.error(
"Cannot find equivalent atoms.\n"
"Multiple atomic species were marked as being equivalent.\n"
"Species marked in this way: {}".format(list(sset))
)
#end if
species_by_specnum[sn] = list(sset)[0]
#end for
# give each unique species a unique label
labels_by_specnum = obj()
species_list = list(species_by_specnum.values())
species_set = set(species_list)
species_counts = obj()
for s in species_set:
species_counts[s] = species_list.count(s)
#end for
spec_counts = obj()
for sn,s in species_by_specnum.items():
if species_counts[s]==1:
labels_by_specnum[sn] = s
else:
if s not in spec_counts:
spec_counts[s] = 1
else:
spec_counts[s] += 1
#end if
labels_by_specnum[sn] = s + str(spec_counts[s])
#end if
#end for
# find indices for each unique species
equiv_indices = obj()
for s in labels_by_specnum.values():
equiv_indices[s] = list()
#end for
for i,sn in enumerate(ds.equivalent_atoms):
equiv_indices[labels_by_specnum[sn]].append(i)
#end for
for s,indices in equiv_indices.items():
equiv_indices[s] = np.array(indices,dtype=int)
#end for
return equiv_indices
#end def equivalent_atoms
# operations to support restricted cases for RMG code
# supported rmg lattices
rmg_lattices = obj(
orthorhombic_P = 'Orthorhombic Primitive',
tetragonal_P = 'Tetragonal Primitive',
hexagonal_P = 'Hexagonal Primitive',
cubic_P = 'Cubic Primitive',
cubic_I = 'Cubic Body Centered',
cubic_F = 'Cubic Face Centered',
)
[docs]
def rmg_lattice(self,allow_tile=False,all_results=False,ret_bravais=False,exit=False,warn=False):
# output variables
rmg_lattice = None
tmatrix = None
rmg_lattices = Structure.rmg_lattices
# represent current bravais lattice
s = Structure(
units = str(self.units),
axes = self.axes.copy(),
elem = ['H'],
pos = [[0,0,0]],
)
# get standard primitive cell based on bravais lattice
sp = s.primitive()
# get current bravais lattice name
d = sp.symmetry_data()
bv = sp.bravais_lattice_name(d)
bvp = bv.rsplit('_',1)[0]+'_P'
if bv in rmg_lattices:
rmg_lattice = bv
#end if
# attempt to get a valid cell by tiling if current one is not valid
if rmg_lattice is None and bvp in rmg_lattices and allow_tile:
spt = Structure(
units = self.units,
axes = d.std_lattice,
)
tmatrix,valid_by_tiling = spt.tilematrix(sp,status=True)
if not valid_by_tiling:
tmatrix = None
else:
# apply tiling matrix
st = sp.copy().tile(tmatrix)
# update lattice type
rmg_lattice,bv = st.rmg_lattice(ret_bravais=True)
#end if
#end if
if rmg_lattice is None and (exit or warn):
msg = (
"Bravais lattice is not supported by the RMG code.\n"
"Cell bravais lattice: {}\n"
"Lattices supported by RMG: {}".format(
bv, list(sorted(rmg_lattices.keys()))
)
)
if exit:
self.error(msg)
elif warn:
self.warn(msg)
#end if
#end if
if all_results:
return rmg_lattice,tmatrix,s,sp,bv
elif allow_tile:
return rmg_lattice,tmatrix
elif ret_bravais:
return rmg_lattice,bv
else:
return rmg_lattice
#end if
#end def rmg_lattice
#end if
#end def rmg_transform
#end class Structure
Structure.set_operations()
#======================#
# SeeK-path functions #
#======================#
# installation instructions for seekpath interface
#
# installation of seekpath
# pip install seekpath
try:
import seekpath
from seekpath import get_explicit_k_path
version = seekpath.__version__
try:
version = [int(i) for i in version.split('.')]
if len(version) < 3:
raise ValueError
#end if
except ValueError:
raise ValueError("Unable to parse version number")
#end try
if tuple(version) < (1, 8, 3):
raise ValueError("Invalid seekpath version, need >= 1.8.4")
#end if
del version
del seekpath
except:
get_explicit_k_path = unavailable('seekpath','get_explicit_k_path')
#end try
def _getseekpath(
structure = None,
with_time_reversal = False,
recipe = 'hpkot',
reference_distance = 0.025,
threshold = 1E-7,
symprec = 1E-5,
angle_tolerance = 1.0,
):
if not isinstance(structure, Structure):
raise TypeError(
"structure is not of type Structure\n"
"type received: {0}".format(structure.__class__.__name__)
)
#end if
if structure.has_folded():
structure = structure.folded_structure
#end if
structure = structure.copy()
if structure.units != 'A':
structure.change_units('A')
#end if
axes = structure.axes
unit_pos = structure.get_scaled_positions()
atomic_num = structure.get_atomic_numbers()
cell = (axes,unit_pos,atomic_num)
return get_explicit_k_path(cell)
#end def _getseekpath
[docs]
def get_conventional_cell(
structure = None,
symprec = 1E-5,
angle_tolerance = 1.0,
seekpathout = None,
):
if seekpathout is None:
seekpathout = _getseekpath(structure=structure, symprec = symprec, angle_tolerance=angle_tolerance)
#end if
try:
axes = seekpathout['conv_lattice']
enumbers = seekpathout['conv_types']
posd = seekpathout['conv_positions']
volfac = seekpathout['volume_original_wrt_conv']
except: # dict-like access above deprecated at some point
axes = seekpathout.conv_lattice
enumbers = seekpathout.conv_types
posd = seekpathout.conv_positions
volfac = seekpathout.volume_original_wrt_conv
bcharge = structure.background_charge*volfac
pos = dot(posd,axes)
sout = structure.copy()
elem = np.array([Elements(i).symbol for i in enumbers], dtype=str)
if abs(bcharge-int(bcharge)) > 1E-6:
raise ValueError("Invalid background charge for conventional structure")
#end if
return {'structure': Structure(axes=axes, elem=elem, pos=pos, background_charge = bcharge, units='A')}
#end def get_conventional_cell
[docs]
def get_primitive_cell(
structure = None,
symprec = 1E-5,
angle_tolerance = 1.0,
seekpathout = None,
):
if seekpathout is None:
seekpathout = _getseekpath(structure = structure, symprec = symprec, angle_tolerance=angle_tolerance)
#end if
axes = seekpathout['primitive_lattice']
enumbers = seekpathout['primitive_types']
posd = seekpathout['primitive_positions']
volfac = seekpathout['volume_original_wrt_prim']
bcharge = structure.background_charge*volfac
pos = dot(posd,axes)
sout = structure.copy()
elem = np.array([Elements(i).symbol for i in enumbers], dtype=str)
return {'structure' : Structure(axes=axes, elem=elem, pos=pos, background_charge=bcharge, units='A'),
'T' : seekpathout['primitive_transformation_matrix']}
#end def get_primitive_cell
[docs]
def get_kpath(
structure = None,
check_standard = True,
with_time_reversal = False,
recipe = 'hpkot',
reference_distance = 0.025,
threshold = 1E-7,
symprec = 1E-5,
angle_tolerance = 1.0,
seekpathout = None,
):
if seekpathout is None:
seekpathout = _getseekpath(structure=structure, symprec = symprec, angle_tolerance=angle_tolerance,
recipe=recipe, reference_distance=reference_distance, with_time_reversal=with_time_reversal)
#end if
if check_standard:
structure = structure.copy()
structure.change_units('A')
axes = structure.axes
primlat = seekpathout['primitive_lattice']
if not np.isclose(primlat, axes).all():
#print primlat, axes
Structure.class_error(
'Input lattice is not the conventional lattice. If you like otherwise, set check_standard=False.'
)
#end if
#end if
inverse_A_to_inverse_B = convert(1.0,'A','B')
return {'explicit_kpoints_abs_inv_A' : seekpathout['explicit_kpoints_abs'],
'explicit_kpoints_abs_inv_B' : seekpathout['explicit_kpoints_abs']*inverse_A_to_inverse_B,
'explicit_kpoints_rel' : seekpathout['explicit_kpoints_rel'],
'explicit_kpoints_labels' : seekpathout['explicit_kpoints_labels'],
'path' : seekpathout['path'],
'explicit_path_linearcoords': seekpathout['explicit_kpoints_linearcoord'],
'point_coords' : seekpathout['point_coords']}
#end def get_kpath
[docs]
def get_symmetry(
structure = None,
symprec = 1E-5,
angle_tolerance = 1.0,
seekpathout = None,
):
if seekpathout is None:
seekpathout = _getseekpath(structure = structure, symprec = symprec, angle_tolerance=angle_tolerance)
#end if
sgint = seekpathout['spacegroup_international']
bravais = seekpathout['bravais_lattice']
invsym = seekpathout['has_inversion_symmetry']
sgnum = seekpathout['spacegroup_number']
return {'sgint': sgint, 'bravais': bravais, 'inv_sym_exists': invsym, 'sgnum': sgnum}
#end def get_symmetry
[docs]
def get_structure_with_bands(
cell = 0,
structure = None,
with_time_reversal = False,
reference_distance = 0.025,
threshold = 1E-7,
symprec = 1E-5,
angle_tolerance = 1.0,
):
if cell == 0:
''' Use input structure '''
struct_band = structure.copy()
elif cell == 1:
''' Use conventional structure '''
struct_band = get_conventional_cell(structure=structure, symprec=symprec, angle_tolerance=angle_tolerance)['structure']
elif cell == 2:
''' Use primitive structure '''
struct_band = get_primitive_cell(structure=structure, symprec=symprec, angle_tolerance=angle_tolerance)['structure']
else:
Structure.class_error('Invalid cell type')
#end if
kpath = get_kpath(structure=struct_band, check_standard=False, with_time_reversal=with_time_reversal)
return Structure(axes = struct_band.axes,
elem = struct_band.elem,
pos = struct_band.pos,
background_charge = struct_band.background_charge,
kpoints = kpath['explicit_kpoints_rel'],
units = 'A')
#end def get_structure_with_bands
# test needed
[docs]
def get_band_tiling(
structure = None,
check_standard = True,
use_ktol = True,
kpoints_label = None,
kpoints_rel = None,
max_volfac = 20,
min_volfac = 0,
target_volfac = None,
):
def cube_deviation(axes):
a = axes
volume = np.abs(dot(cross(axes[0,:], axes[1,:]), axes[2,:]))
dc = volume**(1./3)*sqrt(2.)
d1 = np.abs(norm(a[0]+a[1])-dc)
d2 = np.abs(norm(a[1]+a[2])-dc)
d3 = np.abs(norm(a[2]+a[0])-dc)
d4 = np.abs(norm(a[0]-a[1])-dc)
d5 = np.abs(norm(a[1]-a[2])-dc)
d6 = np.abs(norm(a[2]-a[0])-dc)
return (d1+d2+d3+d4+d5+d6)/(6*dc)
#end def cube_deviation
def cuboid_with_int_edges(vol):
# Given a volume, return the cuboids which have integer edges
divisors = []
edges = []
if isinstance(vol, int):
i = 1
while i<=vol:
if vol%i==0:
divisors.append(i)
#end if
i+=1
#end while
for i in divisors:
for j in divisors:
for k in divisors:
if i*j*k == vol:
edges.append([i,j,k])
#end if
#end for
#end for
#end for
else:
self.error('Volume multiplier must be integer')
#end if
return edges
#end def cuboid_with_int_edges
def alphas_on_grid(alphas, divs):
new_alphas = []
for alpha in alphas:
abs_alpha = np.abs(alpha)
sign_alpha = np.sign(alpha)
new_alpha = np.round(abs_alpha*divs)*1./divs*sign_alpha
new_alphas.append(new_alpha)
#end for
return new_alphas
#end def alphas_on_grid
def find_alphas(structure, kpoints_label, kpoints_rel, check_standard):
# Read wavevectors from the input and return the differences between all wavevectors (alphas)
kpath = get_kpath(structure = structure, check_standard = check_standard)
kpath_label = np.array(kpath['explicit_kpoints_labels'])
kpath_rel = kpath['explicit_kpoints_rel']
kpts = dict()
if kpoints_label is None:
kpoints_label = []
if kpoints_rel is None:
Structure.class_error(
"Please define symbolic or crystal coordinates for kpoints. e.g. ['GAMMA', 'K'] or [[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]]"
)
else:
for k in kpoints_rel:
kindex = np.isclose(kpath_rel,k, atol=1e-5).all(1)
if any(kindex):
kpts[kpath_label[kindex][0]] = np.array(k)
else:
Structure.class_error('{0} is not found in the kpath'.format(k))
#end if
#end for
#end if
else:
if kpoints_rel is not None:
Structure.class_error('Both symbolic and crystal k-points are defined.')
else:
kpoints_rel = []
num_kpoints = 0
for k in kpoints_label:
kindex = k == kpath_label
if any(kindex):
if k == '' or k is None:
k = '{0}'.format(num_kpoints)
#end if
kpts[k] = np.array(kpath_rel[kindex][0])
else:
Structure.class_error('{0} is not found in the kpath'.format(k))
#end if
#end for
#end if
#end if
alphas = np.array([x[0] - x[1] for x in itertools.combinations(kpts.values(),2)]) #Combinations of k_1 - k_2 in kpts list
kpt0 = list(kpts.values())[0]
return alphas, kpt0
#end def find_alphas
def find_vars(alphas,min_volfac, max_volfac, target_volfac, use_ktol):
'''
Find variables to generate possible smallest matrices in the Upper Triangular Hermite Normal Form, from PHYSICAL REVIEW B 92, 184301 (2015)
For target or min/max volumes, it returns vol_mul, which is the smallest volume multiplier to be used on the volfac here
'''
if use_ktol:
ktol = 0.25/max_volfac
else:
ktol = 0.0
#end if
if target_volfac is not None:
if min_volfac is None and max_volfac is None:
min_volfac = target_volfac
max_volfac = target_volfac
else:
print("target_volfac and {min_volfac, max_volfac} cannot be defined together!")
exit()
#end if
#end if
cur_volfac = 1e6
vars = []
for mvol in range(max_volfac, 0, -1):
cuboids = cuboid_with_int_edges(mvol)
for c in cuboids:
a_1, a_2, a_3 = c
new_alphas = alphas_on_grid(alphas, c)
rec_grid_voxel = np.array([1./a_1,1./a_2,1./a_3]) # reciprocal grid voxel size
rem = []
for alpha in alphas:
rem.append(np.mod(np.abs(alpha), rec_grid_voxel))
#end for
if np.all(np.isclose(rem,0.0,atol=ktol)):
n1 = a_1
n2 = a_2
n3 = a_3
g12 = np.gcd.reduce([n1,n2])
g13 = np.gcd.reduce([n1,n3])
g23 = np.gcd.reduce([n2,n3])
g123 = np.gcd.reduce([n1, n2, n3])
volfac = n1*n2*n3*g123//(g12*g13*g23)
if volfac < cur_volfac: #min_volfac <= volfac and and volfac <= max_volfac:
vars = [[n1, n2, n3, g12, g13, g23, g123]]
cur_volfac = volfac
elif volfac == cur_volfac:
vars.append([n1, n2, n3, g12, g13, g23, g123])
#end if
#end if
#end for
#end for
if vars == []:
print('Change ktol')
exit()
else:
can_be_found = False
vol_mul = 1
while not can_be_found:
if volfac*vol_mul <= max_volfac and volfac*vol_mul >= min_volfac:
can_be_found = True
elif volfac*vol_mul > max_volfac:
print('Increase max_volfac or target_volfac!')
exit()
else:
vol_mul+=1
#end if
#end while
return vars, vol_mul
#end if
#end def find_vars
def find_mats(mat_vars, alphas):
'''
Given the variables (v), return the list of all upper triangular matrices as in PHYSICAL REVIEW B 92, 184301 (2015)
'''
mats = []
for v in mat_vars:
n1, n2, n3, g12, g13, g23, g123 = v
#New alphas exactly on the voxels thanks to ktol
for p in range(0, g23):
for q in range(0, g12//g123):
for r in range(g13*g23//g123):
temp_mat = [
[g123*n1//(g12*g13), q*g123*n2//(g12*g23), r*g123*n3//(g13*g23)],
[0, n2//g23, p*n3//g23],
[0,0,n3]]
comm = []
div = np.array([n1, n2, n3])
new_alphas = alphas_on_grid(alphas, div)
for new_alpha in new_alphas:
if (np.isclose(np.abs(dot(temp_mat,new_alpha))%1.0, 0, atol = 1e-6)).all():
comm.append(True) # new_alpha is commensurate with tmat
else:
comm.append(False)
#end if
#end for
if all(comm) and temp_mat not in mats: # if all new_alphas are commensurate with tmat
mats.append(temp_mat)
#end if
#end for
#end for
#end for
#end for
return mats
#end def find_mats
def find_cubic_mat(mats, structure, mat_vol_mul):
final_axes = []
final_t = []
final_cubicity = 1e6
mats = np.array(mats)
for m in mats:
axes = structure.axes.copy()
s = structure.copy()
[m_t, r] = optimal_tilematrix(s.tile(m), volfac=mat_vol_mul)
m_axes = dot(dot(m_t, m), axes)
m_cubicity = cube_deviation(m_axes)
if m_cubicity < final_cubicity:
final_axes = m_axes
final_cubicity = m_cubicity
final_t = dot(m_t, m)
#end if
#end for
return final_t.tolist()
#end def find_cubic_mat
def find_shift(final_mat, structure, kpt0):
return None
#end def find_cubic_mat
alphas, kpt0 = find_alphas(structure,kpoints_label,kpoints_rel, check_standard) # Wavevector differences
mat_vars, mat_vol_mul = find_vars(alphas,min_volfac,max_volfac,target_volfac,use_ktol) # Variables to construct upper triangular matrices
mats = find_mats(mat_vars,alphas) # List of upper triangular matrices that are commensurate with alphas
final_mat = find_cubic_mat(mats, structure, mat_vol_mul) # Matrix leading to a lattice with highest cubicity, optimized using elementary operations
shift = find_shift(final_mat, structure, kpts0) # Find the grid shift
o = obj()
o.mat = final_mat
o.shift = shift
o.det = det(final_mat)
return o
#end def get_band_tiling
[docs]
def get_seekpath_full(
structure = None,
seekpathout = None,
conventional = False,
primitive = False,
**kwargs
):
if seekpathout is None:
seekpathout = _getseekpath(structure,**kwargs)
#end if
res = obj(seekpathout)
for k,v in res.items():
if isinstance(v,dict):
res[k] = obj(v)
#end if
#end for
if conventional:
conv = get_conventional_cell(structure,seekpathout=seekpathout)
res.conventional = conv['structure']
#end if
if primitive:
prim = get_primitive_cell(structure,seekpathout=seekpathout)
res.primitive = prim['structure']
res.prim_tmatrix = prim['T']
#end if
return res
#end def get_seekpath_full
skp = obj(
_getseekpath = _getseekpath,
get_conventional_cell = get_conventional_cell,
get_primitive_cell = get_primitive_cell,
get_kpath = get_kpath,
get_symmetry = get_symmetry,
get_structure_with_bands = get_structure_with_bands,
get_band_tiling = get_band_tiling,
)
#==========================#
# end SeeK-path functions #
#==========================#
[docs]
def interpolate_structures(struct1,struct2=None,images=None,min_image=True,recenter=True,match_com=False,repackage=False,chained=False):
if images is None:
Structure.class_error('images must be provided','interpolate_structures')
#end if
# if a list of structures is provided,
# interpolate between pairs in the chain of structures
if isinstance(struct1,(list,tuple)):
structures_in = struct1
structures = []
for n in range(len(structures_in)-1):
struct1 = structures_in[n]
struct2 = structures_in[n+1]
structs = interpolate_structures(struct1,struct2,images,min_image,recenter,match_com,repackage,chained=True)
if n==0:
structures.append(structs[0])
#end if
structures.extend(structs[1:-1])
if n==len(structures_in)-2:
structures.append(structs[-1])
#end if
#end for
return structures
#end if
# handle PhysicalSystem objects indirectly
system1 = None
system2 = None
if not isinstance(struct1,Structure):
system1 = struct1.copy()
system1.remove_folded()
struct1 = system1.structure
#end if
if not isinstance(struct2,Structure):
system2 = struct2.copy()
system2.remove_folded()
struct2 = system2.structure
#end if
# perform the interpolation
structures = struct1.interpolate(struct2,images,min_image,recenter,match_com)
# repackage into physical system objects if requested
if repackage:
if system1 is not None:
system = system1
elif system2 is not None:
system = system2
else:
Structure.class_error('cannot repackage into physical systems since no system object was provided in place of a structure','interpolate_structures')
#end if
systems = []
for s in structures:
ps = system.copy()
ps.structure = s
systems.append(ps)
#end for
result = systems
else:
result = structures
#end if
return result
#end def interpolate_structures
# test needed
[docs]
def structure_animation(filepath,structures,tiling=None):
path,file = os.path.split(filepath)
if not file.endswith('xyz'):
Structure.class_error('only xyz files are supported for now','structure_animation')
#end if
anim = ''
for s in structures:
if tiling is None:
anim += s.write_xyz()
else:
anim += s.tile(tiling).write_xyz()
#end if
#end for
open(filepath,'w').write(anim)
#end def structure_animation
[docs]
class DefectStructure(Structure):
def __init__(self,*args,**kwargs):
if len(args)>0 and isinstance(args[0],Structure):
self.transfer_from(args[0],copy=True)
else:
Structure.__init__(self,*args,**kwargs)
#end if
#end def __init__
[docs]
def defect_from_bond_compression(self,compression_cutoff,bond_eq,neighbors):
bind,bcent,blens = self.bonds(neighbors)
ind = bind[ np.abs(blens/bond_eq - 1.) > compression_cutoff ]
idefect = np.array(list(set(ind.ravel())))
defect = self.carve(idefect)
return defect
#end def defect_from_bond_compression
[docs]
def defect_from_displacement(self,displacement_cutoff,reference):
displacement = self.scalar_displacement(reference)
idefect = displacement > displacement_cutoff
defect = self.carve(idefect)
return defect
#end def defect_from_displacement
[docs]
def compare(self,dist_cutoff,d1,d2=None):
if d2 is None:
d2 = d1
d1 = self
#end if
res = Sobj()
natoms1 = len(d1.pos)
natoms2 = len(d2.pos)
if natoms1<natoms2:
dsmall,dlarge = d1,d2
else:
dsmall,dlarge = d2,d1
#end if
nn = nearest_neighbors(1,dlarge,dsmall)
dist = dsmall.distances(dlarge[nn.ravel()])
dmatch = dist<dist_cutoff
ismall = np.array(list(range(len(dsmall.pos))))
ismall = ismall[dmatch]
ilarge = nn[ismall]
if natoms1<natoms2:
i1,i2 = ismall,ilarge
else:
i2,i1 = ismall,ilarge
#end if
natoms_match = dmatch.sum()
res.all_match = natoms1==natoms2 and natoms1==natoms_match
res.natoms_match = natoms_match
res.imatch1 = i1
res.imatch2 = i2
return res
#end def compare
#end class DefectStructure
[docs]
class Crystal(Structure):
"""Generate a crystal structure.
Attributes
----------
lattice_constants
lattices
centering_types
lattice_centerings
centerings
cell_types
cell_aliases
cell_classes
constants : NDArray of float
The lattice constants (``a``, ``b``, ``c``) for the crystal.
angles : NDArray of float
The angles (``α``, ``β``, ``γ``) for the crystal.
generation_info : obj of str: str
The supplied inputs to the class constructor.
Parameters
----------
lattice : str, optional
cell : str, optional
centering : {"P", "A", "B", "C", "F", "I", "R"}, optional
constants : float or tuple of float, optional
Lattice constants required for the specified lattice. The order
for these is ``(a, b, c, α, β, γ)``. If the specified lattice
does not require some constant, you can omit it, but retain the
overall order of the constants.
atoms : str or tuple of str, optional
The atomic symbol(s) of the atoms in the lattice.
basis : list of lists of floats, optional
A list of vectors that define the atom positions with respect
to the ``basis_vectors``. If there are multiple ``atoms``, this
should have the sample length as ``atoms``.
basis_vectors : ArrayLike of float or {"prim", "conv"}, optional
A set of 3 vectors that define the basis used to transform
``basis``.
cscale : list of float, optional
Scaling values for the provided constants. Must have the same
length as ``constants``.
See Also
--------
Structure :
All remaining parameters are passed to this class's constructor.
See its docstring for more details.
"""
lattice_constants = obj(
triclinic = ['a','b','c','alpha','beta','gamma'],
monoclinic = ['a','b','c','beta'],
orthorhombic = ['a','b','c'],
tetragonal = ['a','c'],
hexagonal = ['a','c'],
cubic = ['a'],
rhombohedral = ['a','alpha']
)
"""Mapping from a lattice type to the required values to create the cell."""
lattices = list(lattice_constants.keys())
"""List of lattice systems."""
centering_types = obj(
primitive = 'P',
base_centered = ('A','B','C'),
face_centered = 'F',
body_centered = 'I',
rhombohedral_centered = 'R'
)
"""Mapping from centering types to their Pearson symbol."""
lattice_centerings = obj(
triclinic = ['P'],
monoclinic = ['P','A','B','C'],
orthorhombic = ['P','C','I','F'],
tetragonal = ['P','I'],
hexagonal = ['P','R'],
cubic = ['P','I','F'],
rhombohedral = ['P']
)
"""Mapping of lattice systems to allowed centering types."""
centerings = obj(
P = [],
A = [[0,.5,.5]],
B = [[.5,0,.5]],
C = [[.5,.5,0]],
F = [[0,.5,.5],[.5,0,.5],[.5,.5,0]],
I = [[.5,.5,.5]],
R = [[2./3, 1./3, 1./3],[1./3, 2./3, 2./3]]
)
cell_types = set(['primitive','conventional'])
"""Types of cells, currently only ``primitive`` and ``conventional``."""
cell_aliases = obj(
prim = 'primitive',
conv = 'conventional'
)
"""Mapping from shortened aliases ``prim`` and ``conv`` to their cell type."""
cell_classes = obj(
sc = 'cubic',
bcc = 'cubic',
fcc = 'cubic',
hex = 'hexagonal'
)
"""Mapping from common lattice names to their lattices."""
for lattice in lattices:
cell_classes[lattice]=lattice
#end for
#helpful websites for structures
# wikipedia.org
# webelements.com
# webmineral.com
# springermaterials.com
known_crystals = {
('diamond','fcc'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 3.57,
units = 'A',
atoms = 'C',
basis = [[0,0,0],[.25,.25,.25]]
),
('diamond','sc'):obj(
lattice = 'cubic',
cell = 'conventional',
centering = 'F',
constants = 3.57,
units = 'A',
atoms = 'C',
basis = [[0,0,0],[.25,.25,.25]]
),
('diamond','prim'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 3.57,
units = 'A',
atoms = 'C',
basis = [[0,0,0],[.25,.25,.25]]
),
('diamond','conv'):obj(
lattice = 'cubic',
cell = 'conventional',
centering = 'F',
constants = 3.57,
units = 'A',
atoms = 'C',
basis = [[0,0,0],[.25,.25,.25]]
),
('wurtzite','prim'):obj(
lattice = 'hexagonal',
cell = 'primitive',
centering = 'P',
constants = (3.35,5.22),
units = 'A',
#atoms = ('Zn','O'),
#basis = [[1./3, 2./3, 3./8],[1./3, 2./3, 0]]
atoms = ('Zn','O','Zn','O'),
basis = [[0,0,5./8],[0,0,0],[2./3,1./3,1./8],[2./3,1./3,1./2]]
),
('ZnO','prim'):obj(
lattice = 'wurtzite',
cell = 'prim',
constants = (3.35,5.22),
units = 'A',
atoms = ('Zn','O','Zn','O')
),
('NaCl','prim'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 5.64,
units = 'A',
atoms = ('Na','Cl'),
basis = [[0,0,0],[.5,0,0]],
basis_vectors = 'conventional'
),
('rocksalt','prim'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 5.64,
units = 'A',
atoms = ('Na','Cl'),
basis = [[0,0,0],[.5,0,0]],
basis_vectors = 'conventional'
),
('copper','prim'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 3.615,
units = 'A',
atoms = 'Cu'
),
('calcium','prim'):obj(
lattice = 'cubic',
cell = 'primitive',
centering = 'F',
constants = 5.588,
units = 'A',
atoms = 'Ca'
),
# http://www.webelements.com/oxygen/crystal_structure.html
# Phys Rev 160 694
('oxygen','prim'):obj(
lattice = 'monoclinic',
cell = 'primitive',
centering = 'C',
constants = (5.403,3.429,5.086,132.53),
units = 'A',
angular_units = 'degrees',
atoms = ('O','O'),
basis = [[0,0,1.15/2],[0,0,-1.15/2]],
basis_vectors = np.identity(3)
),
# http://en.wikipedia.org/wiki/Calcium_oxide
# http://www.springermaterials.com/docs/info/10681719_224.html
('CaO','prim'):obj(
lattice = 'NaCl',
cell = 'prim',
constants = 4.81,
atoms = ('Ca','O')
),
('CaO','conv'):obj(
lattice = 'NaCl',
cell = 'conv',
constants = 4.81,
atoms = ('Ca','O')
),
# http://en.wikipedia.org/wiki/Copper%28II%29_oxide
# http://iopscience.iop.org/0953-8984/3/28/001/
# http://www.webelements.com/compounds/copper/copper_oxide.html
# http://www.webmineral.com/data/Tenorite.shtml
('CuO','prim'):obj(
lattice = 'monoclinic',
cell = 'primitive',
centering = 'C',
constants = (4.683,3.422,5.128,99.54),
units = 'A',
angular_units = 'degrees',
atoms = ('Cu','O','Cu','O'),
basis = [[.25,.25,0],[0,.418,.25],
[.25,.75,.5],[.5,.5-.418,.75]],
basis_vectors = 'conventional'
),
('Ca2CuO3','prim'):obj(# kateryna foyevtsova
lattice = 'orthorhombic',
cell = 'primitive',
centering = 'I',
constants = (3.77,3.25,12.23),
units = 'A',
atoms = ('Cu','O','O','O','Ca','Ca'),
basis = [[ 0, 0, 0 ],
[ .50, 0, 0 ],
[ 0, 0, .16026165],
[ 0, 0, .83973835],
[ 0, 0, .35077678],
[ 0, 0, .64922322]],
basis_vectors = 'conventional'
),
('La2CuO4','prim'):obj( #tetragonal structure
lattice = 'tetragonal',
cell = 'primitive',
centering = 'I',
constants = (3.809,13.169),
units = 'A',
atoms = ('Cu','O','O','O','O','La','La'),
basis = [[ 0, 0, 0],
[ .5, 0, 0],
[ 0, .5, 0],
[ 0, 0, .182],
[ 0, 0, -.182],
[ 0, 0, .362],
[ 0, 0, -.362]]
),
('Cl2Ca2CuO2','prim'):obj(
lattice = 'tetragonal',
cell = 'primitive',
centering = 'I',
constants = (3.869,15.05),
units = 'A',
atoms = ('Cu','O','O','Ca','Ca','Cl','Cl'),
basis = [[ 0, 0, 0 ],
[ .5, 0, 0 ],
[ 0, .5, 0 ],
[ .5, .5, .104],
[ 0, 0, .396],
[ 0, 0, .183],
[ .5, .5, .317]],
basis_vectors = 'conventional'
),
('Cl2Ca2CuO2','afm'):obj(
lattice = 'tetragonal',
cell = 'conventional',
centering = 'P',
axes = [[.5,-.5,0],[.5,.5,0],[0,0,1]],
constants = (2*3.869,15.05),
units = 'A',
atoms = 4*['Cu','O','O','Ca','Ca','Cl','Cl'],
basis = [[ 0, 0, 0 ], #Cu
[ .25, 0, 0 ],
[ 0, .25, 0 ],
[ .25, .25, .104],
[ 0, 0, .396],
[ 0, 0, .183],
[ .25, .25, .317],
[ .25, .25, .5 ], #Cu
[ .5, .25, .5 ],
[ .25, .5, .5 ],
[ .5, .5, .604],
[ .25, .25, .896],
[ .25, .25, .683],
[ .5, .5, .817],
[ .5, 0, 0 ], #Cu2
[ .75, 0, 0 ],
[ .5, .25, 0 ],
[ .75, .25, .104],
[ .5, 0, .396],
[ .5, 0, .183],
[ .75, .25, .317],
[ .75, .25, .5 ], #Cu2
[ 0, .25, .5 ],
[ .75, .5, .5 ],
[ 0, .5, .604],
[ .75, .25, .896],
[ .75, .25, .683],
[ 0, .5, .817]],
basis_vectors = 'conventional'
),
('CuO2_plane','prim'):obj(
lattice = 'tetragonal',
cell = 'primitive',
centering = 'P',
constants = (3.809,13.169),
units = 'A',
atoms = ('Cu','O','O'),
basis = [[ 0, 0, 0],
[ .5, 0, 0],
[ 0, .5, 0]]
),
('graphite_aa','hex'):obj(
axes = [[1./2,-sqrt(3.)/2,0],[1./2,sqrt(3.)/2,0],[0,0,1]],
constants = (2.462,3.525),
units = 'A',
atoms = ('C','C'),
basis = [[0,0,0],[2./3,1./3,0]]
),
('graphite_ab','hex'):obj(
axes = [[1./2,-sqrt(3.)/2,0],[1./2,sqrt(3.)/2,0],[0,0,1]],
constants = (2.462,3.525),
units = 'A',
cscale = (1,2),
atoms = ('C','C','C','C'),
basis = [[0,0,0],[2./3,1./3,0],[0,0,1./2],[1./3,2./3,1./2]]
),
('graphene','prim'):obj(
lattice = 'hexagonal',
cell = 'primitive',
centering = 'P',
constants = (2.462,15.0),
units = 'A',
atoms = ('C','C'),
basis = [[0,0,0],[2./3,1./3,0]]
),
('graphene','rect'):obj(
lattice = 'orthorhombic',
cell = 'conventional',
centering = 'C',
constants = (2.462,sqrt(3.)*2.462,15.0),
units = 'A',
atoms = ('C','C'),
basis = [[0,0,0],[1./2,1./6,0]]
)
}
"""Mapping from material names and their cell types to their crystal information."""
kc_keys = list(known_crystals.keys())
for (name,cell) in kc_keys:
desc = known_crystals[name,cell]
if cell=='prim' and (name,'conv') not in known_crystals:
cdesc = desc.copy()
if cdesc.cell=='primitive':
cdesc.cell = 'conventional'
known_crystals[name,'conv'] = cdesc
elif cdesc.cell=='prim':
cdesc.cell = 'conv'
known_crystals[name,'conv'] = cdesc
#end if
#end if
#end if
del kc_keys
def __init__(self,
lattice = None,
cell = None,
centering = None,
constants = None,
atoms = None,
basis = None,
basis_vectors = None,
tiling = None,
cscale = None,
axes = None,
units = None,
angular_units = 'degrees',
kpoints = None,
kgrid = None,
mag = None,
frozen = None,
kshift = (0,0,0),
permute = None,
operations = None,
elem = None,
pos = None,
use_prim = None,
add_kpath = False,
symm_kgrid = False,
):
if lattice is None and cell is None and atoms is None and units is None:
return
#end if
gi = obj(
lattice = lattice ,
cell = cell ,
centering = centering ,
constants = constants ,
atoms = atoms ,
basis = basis ,
basis_vectors = basis_vectors ,
tiling = tiling ,
cscale = cscale ,
axes = axes ,
units = units ,
angular_units = angular_units ,
frozen = frozen ,
mag = mag ,
kpoints = kpoints ,
kgrid = kgrid ,
kshift = kshift ,
permute = permute ,
operations = operations ,
elem = elem ,
pos = pos ,
use_prim = use_prim ,
add_kpath = add_kpath ,
symm_kgrid = symm_kgrid ,
)
generation_info = gi.copy()
lattice_in = lattice
if isinstance(lattice,str):
lattice=lattice.lower()
#end if
if isinstance(cell,str):
cell=cell.lower()
#end if
known_crystal = False
if (lattice_in,cell) in self.known_crystals:
known_crystal = True
lattice_info = self.known_crystals[lattice_in,cell].copy()
elif (lattice,cell) in self.known_crystals:
known_crystal = True
lattice_info = self.known_crystals[lattice,cell].copy()
#end if
if known_crystal:
while 'lattice' in lattice_info and 'cell' in lattice_info and (lattice_info.lattice,lattice_info.cell) in self.known_crystals:
li_old = lattice_info
lattice_info = self.known_crystals[li_old.lattice,li_old.cell].copy()
del li_old.lattice
del li_old.cell
lattice_info.transfer_from(li_old,copy=False)
#end while
if 'cell' in lattice_info:
cell = lattice_info.cell
elif cell in self.cell_aliases:
cell = self.cell_aliases[cell]
elif cell in self.cell_classes:
lattice = self.cell_classes[cell]
else:
self.error(
"cell shape "+cell+" is not recognized\n"
" the variable cell_classes or cell_aliases must be updated to include "+cell
)
#end if
if 'lattice' in lattice_info:
lattice = lattice_info.lattice
#end if
if 'angular_units' in lattice_info:
angular_units = lattice_info.angular_units
#end if
inputs = obj(
centering = centering,
constants = constants,
atoms = atoms,
basis = basis,
basis_vectors = basis_vectors,
tiling = tiling,
cscale = cscale,
axes = axes,
units = units
)
for var,val in inputs.items():
if val is None and var in lattice_info:
inputs[var] = lattice_info[var]
#end if
#end for
centering,constants,atoms,basis,basis_vectors,tiling,cscale,axes,units=inputs.list('centering','constants','atoms','basis','basis_vectors','tiling','cscale','axes','units')
#end if
if constants is None:
self.error('the variable constants must be provided')
#end if
if atoms is None:
self.error('the variable atoms must be provided')
#end if
if lattice not in self.lattices:
self.error(
"lattice type "+str(lattice)+" is not recognized\n"
" valid lattice types are: "+str(list(self.lattices))
)
#end if
if cell=='conventional':
if centering is None:
self.error(
"centering must be provided for a conventional cell\n"
" options for a "+lattice+" lattice are: "+str(self.lattice_centerings[lattice])
)
elif centering not in self.centerings:
self.error(
"centering type "+str(centering)+" is not recognized\n"
" options for a "+lattice+" lattice are: "+str(self.lattice_centerings[lattice])
)
#end if
#end if
if isinstance(constants,int) or isinstance(constants,float):
constants=[constants]
#end if
if len(constants)!=len(self.lattice_constants[lattice]):
self.error(
"the "+lattice+" lattice depends on the constants "+str(self.lattice_constants[lattice])+"\n"
" you provided "+str(len(constants))+": "+str(constants)
)
#end if
if isinstance(atoms,str):
if basis is not None:
atoms = len(basis)*[atoms]
else:
atoms=[atoms]
#end if
#end if
if basis is None:
if len(atoms)==1:
basis = [(0,0,0)]
else:
self.error(
"must provide as many basis coordinates as basis atoms\n"
" atoms provided: "+str(atoms)+"\n"
" basis provided: "+str(basis)
)
#end if
#end if
if basis_vectors is not None and not isinstance(basis_vectors,str) and len(basis_vectors)!=3:
self.error(
"3 basis vectors must be given, you provided "+str(len(basis))+":\n"
" "+str(basis_vectors)
)
#end if
if tiling is None:
tiling = (1,1,1)
#end if
if cscale is None:
cscale = len(constants)*[1]
#end if
if len(cscale)!=len(constants):
self.error('cscale and constants must be the same length')
#end if
basis = np.array(basis)
tiling = np.array(tiling,dtype=int)
cscale = np.array(cscale)
constants = cscale*np.array(constants)
a,b,c,alpha,beta,gamma = None,None,None,None,None,None
if angular_units=='radians':
pi_1o2 = pi/2
pi_2o3 = 2*pi/3
elif angular_units=='degrees':
pi_1o2 = 90.
pi_2o3 = 120.
else:
self.error(
"angular units must be radians or degrees\n"
" you provided "+str(angular_units)
)
#end if
if lattice=='triclinic':
a,b,c,alpha,beta,gamma = constants
elif lattice=='monoclinic':
a,b,c,beta = constants
alpha = gamma = pi_1o2
elif lattice=='orthorhombic':
a,b,c = constants
alpha=beta=gamma=pi_1o2
elif lattice=='tetragonal':
a,c = constants
b=a
alpha=beta=gamma=pi_1o2
elif lattice=='hexagonal':
a,c = constants
b=a
alpha=beta=pi_1o2
gamma=pi_2o3
elif lattice=='cubic':
a=constants[0]
b=c=a
alpha=beta=gamma=pi_1o2
elif lattice=='rhombohedral':
a,alpha = constants
b=c=a
beta=gamma=alpha
#end if
if angular_units=='degrees':
alpha *= pi/180
beta *= pi/180
gamma *= pi/180
#end if
points = [[0,0,0]]
#get the conventional axes
sa,ca = sin(alpha),cos(alpha)
sb,cb = sin(beta) ,cos(beta)
sg,cg = sin(gamma),cos(gamma)
y = (ca-cg*cb)/sg
a1c = a*np.array([1,0,0])
a2c = b*np.array([cg,sg,0])
a3c = c*np.array([cb,y,sqrt(sb**2-y**2)])
#a1c = np.array([a,0,0])
#a2c = np.array([b*cos(gamma),b*sin(gamma),0])
#a3c = np.array([c*cos(beta),c*cos(alpha)*sin(beta),c*sin(alpha)*sin(beta)])
axes_conv = np.array([a1c,a2c,a3c]).copy()
if axes is None:
if cell not in self.cell_types:
self.error(
"cell must be primitive or conventional\n"
" You provided: "+str(cell)
)
#end if
if cell=='primitive' and centering=='P':
cell='conventional'
#end if
#get the primitive axes
if centering=='P':
a1 = a1c
a2 = a2c
a3 = a3c
elif centering=='A':
a1 = a1c
a2 = (a2c+a3c)/2
a3 = (-a2c+a3c)/2
elif centering=='B':
a1 = (a1c+a3c)/2
a2 = a2c
a3 = (-a1c+a3c)/2
elif centering=='C':
a1 = (a1c-a2c)/2
a2 = (a1c+a2c)/2
a3 = a3c
elif centering=='I':
a1=[ a/2, b/2,-c/2]
a2=[-a/2, b/2, c/2]
a3=[ a/2,-b/2, c/2]
elif centering=='F':
a1=[a/2, b/2, 0]
a2=[ 0, b/2, c/2]
a3=[a/2, 0, c/2]
elif centering=='R':
a1=[ a, 0, 0]
a2=[ a/2, a*sqrt(3.)/2, 0]
a3=[-a/6, a/(2*sqrt(3.)), c/3]
else:
self.error(
"the variable centering must be specified\n"
" valid options are: P,A,B,C,I,F,R"
)
#end if
axes_prim = np.array([a1,a2,a3])
if cell=='primitive':
axes = axes_prim
elif cell=='conventional':
axes = axes_conv
points.extend(self.centerings[centering])
#end if
elif known_crystal:
axes = dot(np.diag([a,b,c]),np.array(axes))
#end if
points = np.array(points,dtype=float)
elem = []
pos = []
if basis_vectors is None:
basis_vectors = axes
elif isinstance(basis_vectors,str):
if basis_vectors=='primitive':
basis_vectors = axes_prim
elif basis_vectors=='conventional':
basis_vectors = axes_conv
#end if
#end if
nbasis = len(atoms)
for point in points:
for i in range(nbasis):
atom = atoms[i]
bpoint = basis[i]
p = dot(point,axes) + dot(bpoint,basis_vectors)
elem.append(atom)
pos.append(p)
#end for
#end for
pos = np.array(pos)
self.set(
constants = np.array([a,b,c]),
angles = np.array([alpha,beta,gamma]),
generation_info = generation_info
)
Structure.__init__(
self,
axes = axes,
scale = a,
elem = elem,
pos = pos,
center = axes.sum(0)/2,
units = units,
frozen = frozen,
mag = mag,
tiling = tiling,
kpoints = kpoints,
kgrid = kgrid,
kshift = kshift,
permute = permute,
rescale = False,
operations = operations,
use_prim = use_prim,
add_kpath = add_kpath,
symm_kgrid = symm_kgrid,
)
#end def __init__
#end class Crystal
# test needed
[docs]
class Jellium(Structure):
prefactors = obj()
prefactors.transfer_from({1:2*pi,2:4*pi,3:4./3*pi})
def __init__(self,charge=None,background_charge=None,cell=None,volume=None,density=None,rs=None,dim=3,
axes=None,kpoints=None,kweights=None,kgrid=None,kshift=None,units=None,tiling=None):
del tiling
if rs is not None:
if dim not in self.prefactors:
self.error(
"only 1,2, or 3 dimensional jellium is currently supported\n"
" you requested one with dimension {0}".format(dim)
)
#end if
density = 1.0/(self.prefactors[dim]*rs**dim)
#end if
if axes is not None:
cell = axes
#end if
if background_charge is not None:
charge = background_charge
#end if
if cell is not None:
cell = np.array(cell)
dim = len(cell)
volume = det(cell)
elif volume is not None:
volume = float(volume)
cell = volume**(1./dim)*np.identity(dim)
#end if
if density is not None:
density = float(density)
if charge is None and volume is not None:
charge = density*volume
elif volume is None and charge is not None:
volume = charge/density
cell = volume**(1./dim)*np.identity(dim)
#end if
#end if
if charge is None or cell is None:
self.error(
"not enough information to form jellium structure\n"
" information provided:\n"
" charge: {0}\n"
" cell: {1}\n"
" volume: {2}\n"
" density: {3}\n"
" rs: {4}\n"
" dim: {5}".format(charge, cell, volume, density, rs, dim)
)
#end if
Structure.__init__(self,background_charge=charge,axes=cell,dim=dim,kpoints=kpoints,kweights=kweights,kgrid=kgrid,kshift=kshift,units=units)
#end def __init__
[docs]
def density(self):
return self.background_charge/self.volume()
#end def density
[docs]
def rs(self):
return 1.0/(self.density()*self.prefactors[self.dim])**(1./self.dim)
#end def rs
[docs]
def tile(self):
self.not_implemented()
#end def tile
#end class Jellium
# test needed
[docs]
def generate_cell(shape,tiling=None,scale=1.,units=None,struct_type=Structure):
if tiling is None:
tiling = (1,1,1)
#end if
axes = Sobj()
axes.sc = 1.*np.array([[ 1,0,0],[0, 1,0],[0,0, 1]])
axes.bcc = .5*np.array([[-1,1,1],[1,-1,1],[1,1,-1]])
axes.fcc = .5*np.array([[ 1,1,0],[1, 0,1],[0,1, 1]])
ax = dot(np.diag(tiling),axes[shape])
center = ax.sum(0)/2
c = Structure(axes=ax,scale=scale,center=center,units=units)
if struct_type!=Structure:
c=c.upcast(struct_type)
#end if
return c
#end def generate_cell
[docs]
def generate_structure(type='crystal',*args,**kwargs):
if type=='crystal':
s = generate_crystal_structure(*args,**kwargs)
elif type=='defect':
s = generate_defect_structure(*args,**kwargs)
elif type=='atom':
s = generate_atom_structure(*args,**kwargs)
elif type=='dimer':
s = generate_dimer_structure(*args,**kwargs)
elif type=='trimer':
s = generate_trimer_structure(*args,**kwargs)
elif type=='jellium':
s = generate_jellium_structure(*args,**kwargs)
elif type=='empty':
s = Structure()
elif type=='basic':
s = Structure(*args,**kwargs)
else:
Structure.class_error(
str(type)+" is not a valid structure type\n"
"options are crystal, defect, atom, dimer, trimer, jellium, empty, or basic"
)
#end if
return s
#end def generate_structure
[docs]
def generate_atom_structure(
atom = None,
units = 'A',
Lbox = None,
skew = 0,
axes = None,
kgrid = (1,1,1),
kshift = (0,0,0),
bconds = tuple('nnn'),
struct_type = Structure
):
"""Create a structure with a single atom in the center of a unit cell.
Parameters
----------
atom : str
The atomic symbol of the atom.
units : str, default="A"
The units of the structure, defaults to Angstroms.
Lbox : int or float, optional
Length of the simulation box. Overrides ``axes``.
skew : int or float less than 1, default=0
...
axes : ArrayLike of float, optional
The unit cell axes.
kgrid : tuple of int, default=(1,1,1)
Number of k-points in each direction. Used to create a
Monkhorst-Pack k-point mesh.
kshift : tuple of int, default=(0,0,0)
Vector to use to translate the k-points in the mesh.
bconds : tuple of str, default=("n","n","n")
Boundary conditions for the resulting structure.
"""
if atom is None:
Structure.class_error('atom must be provided','generate_atom_structure')
#end if
if Lbox is not None:
axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]]
#end if
if axes is None:
s = Structure(elem=[atom],pos=[[0,0,0]],units=units,bconds=bconds)
else:
s = Structure(elem=[atom],pos=[[0,0,0]],axes=axes,kgrid=kgrid,kshift=kshift,bconds=bconds,units=units)
s.center_molecule()
#end if
return s
#end def generate_atom_structure
[docs]
def generate_dimer_structure(
dimer = None,
units = 'A',
separation = None,
Lbox = None,
skew = 0,
axes = None,
kgrid = (1,1,1),
kshift = (0,0,0),
bconds = tuple('nnn'),
struct_type = Structure,
axis = 'x'
):
"""Create a structure with a dimer in the center of a unit cell.
Parameters
----------
dimer : list of str
The atomic symbols of the atoms in the dimer.
units : str, default="A"
The units of the structure, defaults to Angstroms.
separation : int or float
The separation between the atoms in the dimer.
Lbox : int or float, optional
Length of the simulation box. Overrides ``axes``.
skew : int or float less than 1, default=0
...
axes : ArrayLike of float, optional
The unit cell axes.
kgrid : tuple of int, default=(1,1,1)
Number of k-points in each direction. Used to create a
Monkhorst-Pack k-point mesh.
kshift : tuple of int, default=(0,0,0)
Vector to use to translate the k-points in the mesh.
bconds : tuple of str, default=("n","n","n")
Boundary conditions for the resulting structure.
axis : {"x", "y", "z"}, optional
The axis that the dimer is aligned on.
"""
if dimer is None:
Structure.class_error('dimer atoms must be provided to construct dimer','generate_dimer_structure')
#end if
if separation is None:
Structure.class_error('separation must be provided to construct dimer','generate_dimer_structure')
#end if
if Lbox is not None:
axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]]
#end if
if axis=='x':
p2 = [separation,0,0]
elif axis=='y':
p2 = [0,separation,0]
elif axis=='z':
p2 = [0,0,separation]
else:
Structure.class_error(
"dimer orientation axis must be x,y,z\n"
" you provided: {0}".format(axis),
"generate_dimer_structure"
)
#end if
if axes is None:
s = Structure(elem=dimer,pos=[[0,0,0],p2],units=units,bconds=bconds)
else:
s = Structure(elem=dimer,pos=[[0,0,0],p2],axes=axes,kgrid=kgrid,kshift=kshift,units=units,bconds=bconds)
s.center_molecule()
#end if
return s
#end def generate_dimer_structure
[docs]
def generate_trimer_structure(
trimer = None,
units = 'A',
separation = None,
angle = None,
Lbox = None,
skew = 0,
axes = None,
kgrid = (1,1,1),
kshift = (0,0,0),
struct_type = Structure,
axis = 'x',
axis2 = 'y',
angular_units = 'degrees',
plane_rot = None
):
"""Create a structure with a dimer in the center of a unit cell.
Parameters
----------
trimer : list of str
The atomic symbols of the atoms in the trimer.
units : str, default="A"
The units of the structure, defaults to Angstroms.
separation : list of ints or floats
The separation between the atoms in the trimer. The first value
is the distance between atom 1 and atom 2, and the second value
is the distance between atom 1 and atom 3.
angle : int or float
The angle formed by the three atoms in the trimer.
Lbox : int or float, optional
Length of the simulation box. Overrides ``axes``.
skew : int or float less than 1, default=0
...
axes : ArrayLike of float, optional
The unit cell axes.
kgrid : tuple of int, default=(1,1,1)
Number of k-points in each direction. Used to create a
Monkhorst-Pack k-point mesh.
kshift : tuple of int, default=(0,0,0)
Vector to use to translate the k-points in the mesh.
axis : {"x", "y", "z"}, optional
The axis that atom 1 and atom 2 of the trimer is aligned on.
axis2 : {"x", "y", "z"}, optional
The axis that atom 1 and atom 3 of the trimer is aligned on.
angular_units : {"degrees", "rad", "radians"}, optional
The units of the supplied angle.
"""
if trimer is None:
Structure.class_error('trimer atoms must be provided to construct trimer','generate_trimer_structure')
#end if
if separation is None:
Structure.class_error('separation must be provided to construct trimer','generate_trimer_structure')
#end if
if len(separation)!=2:
Structure.class_error(
"two separation distances (atom1-atom2,atom1-atom3) must be provided to construct trimer\n"
"you provided {0} separation distances".format(len(separation)),
'generate_trimer_structure'
)
#end if
if angle is None:
Structure.class_error('angle must be provided to construct trimer','generate_trimer_structure')
#end if
if angular_units=='degrees':
angle *= pi/180
elif not angular_units.startswith('rad'):
Structure.class_error(
"angular units must be degrees or radians\n"
"you provided: {0}".format(angular_units),
'generate_trimer_structure'
)
#end if
if axis==axis2:
Structure.class_error(
"axis and axis2 must be different to define the trimer plane\n"
"you provided {0} for both".format(axis),
'generate_trimer_structure'
)
#end if
if Lbox is not None:
axes = [[Lbox*(1-skew),0,0],[0,Lbox,0],[0,0,Lbox*(1+skew)]]
#end if
p1 = [0,0,0]
if axis=='x':
p2 = [separation[0],0,0]
elif axis=='y':
p2 = [0,separation[0],0]
elif axis=='z':
p2 = [0,0,separation[0]]
else:
Structure.class_error(
"trimer bond1 (atom2-atom1) orientation axis must be x,y,z\n"
" you provided: {0}".format(axis),
'generate_trimer_structure'
)
#end if
r = separation[1]
c = cos(angle)
s = sin(angle)
axpair = axis+axis2
if axpair=='xy':
p3 = [r*c,r*s,0]
elif axpair=='yx':
p3 = [r*s,r*c,0]
elif axpair=='yz':
p3 = [0,r*c,r*s]
elif axpair=='zy':
p3 = [0,r*s,r*c]
elif axpair=='zx':
p3 = [r*s,0,r*c]
elif axpair=='xz':
p3 = [r*c,0,r*s]
else:
Structure.class_error(
"trimer bond2 (atom3-atom1) orientation axis must be x,y,z\n"
" you provided: {0}".format(axis2),
'generate_trimer_structure'
)
#end if
if axes is None:
s = Structure(elem=trimer,pos=[p1,p2,p3],units=units)
else:
s = Structure(elem=trimer,pos=[p1,p2,p3],axes=axes,kgrid=kgrid,kshift=kshift,units=units)
s.center_molecule()
#end if
if plane_rot is not None:
s.rotate_plane(axpair,plane_rot,angular_units)
#end if
return s
#end def generate_trimer_structure
# test needed
[docs]
def generate_jellium_structure(*args,**kwargs):
return Jellium(*args,**kwargs)
#end def generate_jellium_structure
[docs]
def generate_crystal_structure(
lattice = None,
cell = None,
centering = None,
constants = None,
atoms = None,
basis = None,
basis_vectors = None,
tiling = None,
cscale = None,
axes = None,
units = None,
angular_units = 'degrees',
mag = None,
kpoints = None,
kweights = None,
kgrid = None,
kshift = (0,0,0),
permute = None,
operations = None,
struct_type = Crystal,
elem = None,
pos = None,
frozen = None,
posu = None,
elem_pos = None,
folded_elem = None,
folded_pos = None,
folded_units = None,
use_prim = None,
add_kpath = False,
symm_kgrid = False,
#legacy inputs
structure = None,
shape = None,
element = None,
scale = None,
):
"""Generate a crystal structure.
See :py:class:`~.Crystal` and :py:class:`~.Structure` for a
description of the available parameters.
"""
if structure is not None:
lattice = structure
#end if
if shape is not None:
cell = shape
#end if
if element is not None:
atoms = element
#end if
if scale is not None:
constants = scale
#end if
#interface for total manual specification
# this is only here because 'crystal' is default and must handle other cases
s = None
has_elem_and_pos = elem is not None and (pos is not None or posu is not None)
has_elem_and_pos |= elem_pos is not None
if has_elem_and_pos:
s = Structure(
axes = axes,
elem = elem,
pos = pos,
units = units,
mag = mag,
frozen = frozen,
tiling = tiling,
kpoints = kpoints,
kgrid = kgrid,
kshift = kshift,
permute = permute,
rescale = False,
operations = operations,
posu = posu,
elem_pos = elem_pos,
use_prim = use_prim,
add_kpath = add_kpath,
symm_kgrid = symm_kgrid,
)
elif isinstance(structure,Structure):
s = structure
if use_prim is not None and use_prim is not False:
s.become_primitive(source=use_prim,add_kpath=add_kpath)
#end if
if tiling is not None:
s = s.tile(tiling)
#end if
if kpoints is not None:
s.add_kpoints(kpoints,kweights)
#end if
if kgrid is not None:
if not symm_kgrid:
s.add_kmesh(kgrid,kshift)
else:
s.add_symmetrized_kmesh(kgrid,kshift)
#end if
#end if
#end if
if s is not None:
# add point group folded molecular system if present
if folded_elem is not None and folded_pos is not None:
if folded_units is None:
folded_units = units
#end if
fs = Structure(
elem = folded_elem,
pos = folded_pos,
units = folded_units,
rescale = False,
)
s.set_folded(fs)
#end if
return s
#end if
s=Crystal(
lattice = lattice ,
cell = cell ,
centering = centering ,
constants = constants ,
atoms = atoms ,
basis = basis ,
basis_vectors = basis_vectors ,
tiling = tiling ,
cscale = cscale ,
axes = axes ,
units = units ,
angular_units = angular_units ,
frozen = frozen ,
mag = mag ,
kpoints = kpoints ,
kgrid = kgrid ,
kshift = kshift ,
permute = permute ,
operations = operations ,
elem = elem ,
pos = pos ,
use_prim = use_prim ,
add_kpath = add_kpath ,
symm_kgrid = symm_kgrid ,
)
if struct_type!=Crystal:
s=s.upcast(struct_type)
#end if
return s
#end def generate_crystal_structure
defects = obj(
diamond = obj(
H = obj(
pristine = [[0,0,0]],
defect = [[0,0,0],[.625,.375,.375]]
),
T = obj(
pristine = [[0,0,0]],
defect = [[0,0,0],[.5,.5,.5]]
),
X = obj(
pristine = [[.25,.25,.25]],
defect = [[.39,.11,.15],[.11,.39,.15]]
),
FFC = obj(
#pristine = [[ 0, 0, 0],[.25,.25,.25]],
#defect = [[.151,.151,-.08],[.10,.10,.33]]
pristine = [[ 0, 0, 0],[.25 ,.25 ,.25 ],[.5 ,.5 , 0],[.75 ,.75 ,.25 ],[1.5 ,1.5 ,0 ],[1.75 ,1.75 ,.25 ]],
defect = [[.151,.151,-.081],[.099,.099,.331],[.473,.473,-.059],[.722,.722,.230],[1.528,1.528,.020],[1.777,1.777,.309]]
)
)
)
[docs]
def generate_defect_structure(defect,structure,shape=None,element=None,
tiling=None,scale=1.,kgrid=None,kshift=(0,0,0),
units=None,struct_type=DefectStructure):
if structure in defects:
dstruct = defects[structure]
else:
DefectStructure.class_error('defects for '+structure+' structure have not yet been implemented')
#end if
if defect in dstruct:
drep = dstruct[defect]
else:
DefectStructure.class_error(defect+' defect not found for '+structure+' structure')
#end if
ds = generate_crystal_structure(
structure = structure,
shape = shape,
element = element,
tiling = tiling,
scale = 1.0,
kgrid = kgrid,
kshift = kshift,
units = units,
struct_type = struct_type
)
ds.replace(drep.pristine,pos=drep.defect)
ds.rescale(scale)
return ds
#end def generate_defect_structure
[docs]
def read_structure(filepath,elem=None,format=None):
s = generate_structure('empty')
s.read(filepath,elem=elem,format=format)
return s
#end def read_structure
if __name__=='__main__':
large = generate_structure(
type = 'crystal',
structure = 'diamond',
cell = 'fcc',
atoms = 'Ge',
constants = 5.639,
units = 'A',
tiling = (2,2,2),
kgrid = (1,1,1),
kshift = (0,0,0),
)
small = large.folded_structure
print(small.kpoints_unit())
prim = read_structure('scf.struct.xsf')
prim = get_primitive_cell(structure=prim)['structure']
tiling = get_band_tiling(structure=prim, kpoints_label = ['L', 'F'], min_volfac=6, max_volfac = 6)
print(tiling)
#end if
