You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
314 lines
12 KiB
314 lines
12 KiB
from math import sqrt |
|
|
|
from cmmde_atoms import Atoms |
|
from cmmde_symbols import string2symbols |
|
from cmmde_data import reference_states, atomic_numbers, chemical_symbols |
|
from cmmde_utils import plural |
|
|
|
|
|
def incompatible_cell(*, want, have): |
|
return RuntimeError('Cannot create {} cell for {} structure' |
|
.format(want, have)) |
|
|
|
|
|
def bulk(name, crystalstructure=None, a=None, b=None, c=None, *, alpha=None, |
|
covera=None, u=None, orthorhombic=False, cubic=False, |
|
basis=None): |
|
"""Creating bulk systems. |
|
|
|
Crystal structure and lattice constant(s) will be guessed if not |
|
provided. |
|
|
|
name: str |
|
Chemical symbol or symbols as in 'MgO' or 'NaCl'. |
|
crystalstructure: str |
|
Must be one of sc, fcc, bcc, tetragonal, bct, hcp, rhombohedral, |
|
orthorhombic, mlc, diamond, zincblende, rocksalt, cesiumchloride, |
|
fluorite or wurtzite. |
|
a: float |
|
Lattice constant. |
|
b: float |
|
Lattice constant. If only a and b is given, b will be interpreted |
|
as c instead. |
|
c: float |
|
Lattice constant. |
|
alpha: float |
|
Angle in degrees for rhombohedral lattice. |
|
covera: float |
|
c/a ratio used for hcp. Default is ideal ratio: sqrt(8/3). |
|
u: float |
|
Internal coordinate for Wurtzite structure. |
|
orthorhombic: bool |
|
Construct orthorhombic unit cell instead of primitive cell |
|
which is the default. |
|
cubic: bool |
|
Construct cubic unit cell if possible. |
|
""" |
|
|
|
if c is None and b is not None: |
|
# If user passes (a, b) positionally, we want it as (a, c) instead: |
|
c, b = b, c |
|
|
|
if covera is not None and c is not None: |
|
raise ValueError("Don't specify both c and c/a!") |
|
|
|
xref = None |
|
ref = {} |
|
|
|
if name in chemical_symbols: |
|
Z = atomic_numbers[name] |
|
ref = reference_states[Z] |
|
if ref is not None: |
|
xref = ref['symmetry'] |
|
|
|
# If user did not specify crystal structure, and no basis |
|
# is given, and the reference state says we need one, but |
|
# does not have one, then we can't proceed. |
|
if (crystalstructure is None and basis is None |
|
and 'basis' in ref and ref['basis'] is None): |
|
# XXX This is getting much too complicated, we need to split |
|
# this function up. A lot. |
|
raise RuntimeError('This structure requires an atomic basis') |
|
|
|
if ref is None: |
|
ref = {} # easier to 'get' things from empty dictionary than None |
|
|
|
if xref == 'cubic': |
|
# P and Mn are listed as 'cubic' but the lattice constants |
|
# are 7 and 9. They must be something other than simple cubic |
|
# then. We used to just return the cubic one but that must |
|
# have been wrong somehow. --askhl |
|
raise RuntimeError('Only simple cubic ("sc") supported') |
|
|
|
# Mapping of name to number of atoms in primitive cell. |
|
structures = {'sc': 1, 'fcc': 1, 'bcc': 1, |
|
'tetragonal': 1, |
|
'bct': 1, |
|
'hcp': 1, |
|
'rhombohedral': 1, |
|
'orthorhombic': 1, |
|
'mcl': 1, |
|
'diamond': 1, |
|
'zincblende': 2, 'rocksalt': 2, 'cesiumchloride': 2, |
|
'fluorite': 3, 'wurtzite': 2} |
|
|
|
if crystalstructure is None: |
|
crystalstructure = xref |
|
if crystalstructure not in structures: |
|
raise ValueError('No suitable reference data for bulk {}.' |
|
' Reference data: {}' |
|
.format(name, ref)) |
|
|
|
if crystalstructure not in structures: |
|
raise ValueError('Unknown structure: {}.' |
|
.format(crystalstructure)) |
|
|
|
# Check name: |
|
natoms = len(string2symbols(name)) |
|
natoms0 = structures[crystalstructure] |
|
if natoms != natoms0: |
|
raise ValueError('Please specify {} for {} and not {}' |
|
.format(plural(natoms0, 'atom'), |
|
crystalstructure, natoms)) |
|
|
|
if alpha is None: |
|
alpha = ref.get('alpha') |
|
|
|
if a is None: |
|
if xref != crystalstructure: |
|
raise ValueError('You need to specify the lattice constant.') |
|
try: |
|
a = ref['a'] |
|
except KeyError: |
|
raise KeyError('No reference lattice parameter "a" for "{}"' |
|
.format(name)) |
|
|
|
if b is None: |
|
bovera = ref.get('b/a') |
|
if bovera is not None and a is not None: |
|
b = bovera * a |
|
|
|
if crystalstructure in ['hcp', 'wurtzite']: |
|
if cubic: |
|
raise incompatible_cell(want='cubic', have=crystalstructure) |
|
|
|
if c is not None: |
|
covera = c / a |
|
elif covera is None: |
|
if xref == crystalstructure: |
|
covera = ref['c/a'] |
|
else: |
|
covera = sqrt(8 / 3) |
|
|
|
if covera is None: |
|
covera = ref.get('c/a') |
|
if c is None and covera is not None: |
|
c = covera * a |
|
|
|
if orthorhombic and crystalstructure not in ['sc', 'tetragonal', |
|
'orthorhombic']: |
|
return _orthorhombic_bulk(name, crystalstructure, a, covera, u) |
|
|
|
if cubic and crystalstructure in ['bcc', 'cesiumchloride']: |
|
return _orthorhombic_bulk(name, crystalstructure, a, covera) |
|
|
|
if cubic and crystalstructure != 'sc': |
|
return _cubic_bulk(name, crystalstructure, a) |
|
|
|
if crystalstructure == 'sc': |
|
atoms = Atoms(name, cell=(a, a, a), pbc=True) |
|
elif crystalstructure == 'fcc': |
|
b = a / 2 |
|
atoms = Atoms(name, cell=[(0, b, b), (b, 0, b), (b, b, 0)], pbc=True) |
|
elif crystalstructure == 'bcc': |
|
b = a / 2 |
|
atoms = Atoms(name, cell=[(-b, b, b), (b, -b, b), (b, b, -b)], |
|
pbc=True) |
|
elif crystalstructure == 'hcp': |
|
atoms = Atoms(2 * name, |
|
scaled_positions=[(0, 0, 0), |
|
(1 / 3, 2 / 3, 0.5)], |
|
cell=[(a, 0, 0), |
|
(-a / 2, a * sqrt(3) / 2, 0), |
|
(0, 0, covera * a)], |
|
pbc=True) |
|
elif crystalstructure == 'diamond': |
|
atoms = bulk(2 * name, 'zincblende', a) |
|
elif crystalstructure == 'zincblende': |
|
s1, s2 = string2symbols(name) |
|
atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a) |
|
atoms.positions[1] += a / 4 |
|
elif crystalstructure == 'rocksalt': |
|
s1, s2 = string2symbols(name) |
|
atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a) |
|
atoms.positions[1, 0] += a / 2 |
|
elif crystalstructure == 'cesiumchloride': |
|
s1, s2 = string2symbols(name) |
|
atoms = bulk(s1, 'sc', a) + bulk(s2, 'sc', a) |
|
atoms.positions[1, :] += a / 2 |
|
elif crystalstructure == 'fluorite': |
|
s1, s2, s3 = string2symbols(name) |
|
atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a) + bulk(s3, 'fcc', a) |
|
atoms.positions[1, :] += a / 4 |
|
atoms.positions[2, :] += a * 3 / 4 |
|
elif crystalstructure == 'wurtzite': |
|
u = u or 0.25 + 1 / 3 / covera**2 |
|
atoms = Atoms(2 * name, |
|
scaled_positions=[(0, 0, 0), |
|
(1 / 3, 2 / 3, 0.5 - u), |
|
(1 / 3, 2 / 3, 0.5), |
|
(0, 0, 1 - u)], |
|
cell=[(a, 0, 0), |
|
(-a / 2, a * sqrt(3) / 2, 0), |
|
(0, 0, a * covera)], |
|
pbc=True) |
|
elif crystalstructure == 'bct': |
|
from ase.lattice import BCT |
|
if basis is None: |
|
basis = ref.get('basis') |
|
if basis is not None: |
|
natoms = len(basis) |
|
lat = BCT(a=a, c=c) |
|
atoms = Atoms([name] * natoms, cell=lat.tocell(), pbc=True, |
|
scaled_positions=basis) |
|
elif crystalstructure == 'rhombohedral': |
|
atoms = _build_rhl(name, a, alpha, basis) |
|
elif crystalstructure == 'orthorhombic': |
|
atoms = Atoms(name, cell=[a, b, c], pbc=True) |
|
else: |
|
raise ValueError('Unknown crystal structure: ' + crystalstructure) |
|
|
|
if orthorhombic: |
|
assert atoms.cell.orthorhombic |
|
if cubic: |
|
assert abs(atoms.cell.angles() - 90).all() < 1e-10 |
|
return atoms |
|
|
|
|
|
def _build_rhl(name, a, alpha, basis): |
|
from ase.lattice import RHL |
|
lat = RHL(a, alpha) |
|
cell = lat.tocell() |
|
if basis is None: |
|
# RHL: Given by A&M as scaled coordinates "x" of cell.sum(0): |
|
basis_x = reference_states[atomic_numbers[name]]['basis_x'] |
|
basis = basis_x[:, None].repeat(3, axis=1) |
|
natoms = len(basis) |
|
return Atoms([name] * natoms, cell=cell, scaled_positions=basis, pbc=True) |
|
|
|
|
|
def _orthorhombic_bulk(name, crystalstructure, a, covera=None, u=None): |
|
if crystalstructure == 'fcc': |
|
b = a / sqrt(2) |
|
atoms = Atoms(2 * name, cell=(b, b, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)]) |
|
elif crystalstructure == 'bcc': |
|
atoms = Atoms(2 * name, cell=(a, a, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)]) |
|
elif crystalstructure == 'hcp': |
|
atoms = Atoms(4 * name, |
|
cell=(a, a * sqrt(3), covera * a), |
|
scaled_positions=[(0, 0, 0), |
|
(0.5, 0.5, 0), |
|
(0.5, 1 / 6, 0.5), |
|
(0, 2 / 3, 0.5)], |
|
pbc=True) |
|
elif crystalstructure == 'diamond': |
|
atoms = _orthorhombic_bulk(2 * name, 'zincblende', a) |
|
elif crystalstructure == 'zincblende': |
|
s1, s2 = string2symbols(name) |
|
b = a / sqrt(2) |
|
atoms = Atoms(2 * name, cell=(b, b, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0, 0.25), |
|
(0.5, 0.5, 0.5), (0, 0.5, 0.75)]) |
|
elif crystalstructure == 'rocksalt': |
|
s1, s2 = string2symbols(name) |
|
b = a / sqrt(2) |
|
atoms = Atoms(2 * name, cell=(b, b, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0.5, 0), |
|
(0.5, 0.5, 0.5), (0, 0, 0.5)]) |
|
elif crystalstructure == 'cesiumchloride': |
|
atoms = Atoms(name, cell=(a, a, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)]) |
|
elif crystalstructure == 'wurtzite': |
|
u = u or 0.25 + 1 / 3 / covera**2 |
|
atoms = Atoms(4 * name, |
|
cell=(a, a * 3**0.5, covera * a), |
|
scaled_positions=[(0, 0, 0), |
|
(0, 1 / 3, 0.5 - u), |
|
(0, 1 / 3, 0.5), |
|
(0, 0, 1 - u), |
|
(0.5, 0.5, 0), |
|
(0.5, 5 / 6, 0.5 - u), |
|
(0.5, 5 / 6, 0.5), |
|
(0.5, 0.5, 1 - u)], |
|
pbc=True) |
|
else: |
|
raise incompatible_cell(want='orthorhombic', have=crystalstructure) |
|
|
|
return atoms |
|
|
|
|
|
def _cubic_bulk(name, crystalstructure, a): |
|
if crystalstructure == 'fcc': |
|
atoms = Atoms(4 * name, cell=(a, a, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0, 0.5, 0.5), |
|
(0.5, 0, 0.5), (0.5, 0.5, 0)]) |
|
elif crystalstructure == 'diamond': |
|
atoms = _cubic_bulk(2 * name, 'zincblende', a) |
|
elif crystalstructure == 'zincblende': |
|
atoms = Atoms(4 * name, cell=(a, a, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.25, 0.25, 0.25), |
|
(0, 0.5, 0.5), (0.25, 0.75, 0.75), |
|
(0.5, 0, 0.5), (0.75, 0.25, 0.75), |
|
(0.5, 0.5, 0), (0.75, 0.75, 0.25)]) |
|
elif crystalstructure == 'rocksalt': |
|
atoms = Atoms(4 * name, cell=(a, a, a), pbc=True, |
|
scaled_positions=[(0, 0, 0), (0.5, 0, 0), |
|
(0, 0.5, 0.5), (0.5, 0.5, 0.5), |
|
(0.5, 0, 0.5), (0, 0, 0.5), |
|
(0.5, 0.5, 0), (0, 0.5, 0)]) |
|
else: |
|
raise incompatible_cell(want='cubic', have=crystalstructure) |
|
|
|
return atoms
|
|
|