# Source code for ase.ga.bulk_utilities

import numpy as np
from ase.geometry.cell import cell_to_cellpar

def get_cell_angles_lengths(cell):
'''
Returns cell vectors lengths (a,b,c) as well as different
angles (alpha, beta, gamma, phi, chi, psi) (in radians).
'''
cellpar = cell_to_cellpar(cell)
cellpar[3:] *= np.pi / 180  # convert angles to radians
parnames = ['a', 'b', 'c', 'alpha', 'beta', 'gamma']
values = {n: p for n, p in zip(parnames, cellpar)}

volume = abs(np.linalg.det(cell))
for i, param in enumerate(['phi', 'chi', 'psi']):
ab = np.linalg.norm(
np.cross(cell[(i + 1) % 3, :], cell[(i + 2) % 3, :]))
c = np.linalg.norm(cell[i, :])
s = np.abs(volume / (ab * c))
if 1 + 1e-6 > s > 1:
s = 1.
values[param] = np.arcsin(s)

return values

[docs]class CellBounds:
'''
Class for defining as well as checking limits on
cell vector lengths and angles

Parameters:

bounds: dict
Any of the following keywords can be used, in
conjunction with a [low, high] list determining
the lower and upper bounds:

a, b, c:
Minimal and maximal lengths (in Angstrom)
for the 1st, 2nd and 3rd lattice vectors.

alpha, beta, gamma:
Minimal and maximal values (in degrees)
for the angles between the lattice vectors.

phi, chi, psi:
Minimal and maximal values (in degrees)
for the angles between each lattice vector
and the plane defined by the other two vectors.

Example:

>>> from ase.ga.bulk_utilities import CellBounds
>>> CellBounds(bounds={'phi': [20, 160],
...                    'chi': [60, 120],
...                    'psi': [20, 160],
...                    'a': [2, 20], 'b': [2, 20], 'c': [2, 20]})
'''

def __init__(self, bounds={}):
self.bounds = {'alpha': [0, np.pi], 'beta': [0, np.pi],
'gamma': [0, np.pi], 'phi': [0, np.pi],
'chi': [0, np.pi], 'psi': [0, np.pi],
'a': [0, 1e6], 'b': [0, 1e6], 'c': [0, 1e6]}

for param, bound in bounds.items():
if param not in ['a', 'b', 'c']:
# Convert angle from degree to radians
bound = [b * np.pi / 180. for b in bound]
self.bounds[param] = bound

def is_within_bounds(self, cell):
values = get_cell_angles_lengths(cell)
verdict = True
for param, bound in self.bounds.items():
if not (bound[0] <= values[param] <= bound[1]):
verdict = False
return verdict

def get_rotation_matrix(u, t):
'''
Returns the transformation matrix for rotation over an angle t
along an axis with direction u.
'''
ux, uy, uz = u
cost, sint = np.cos(t), np.sin(t)
rotmat = np.array([[(ux**2) * (1 - cost) + cost,
ux * uy * (1 - cost) - uz * sint,
ux * uz * (1 - cost) + uy * sint],
[ux * uy * (1 - cost) + uz * sint,
(uy**2) * (1 - cost) + cost,
uy * uz * (1 - cost) - ux * sint],
[ux * uz * (1 - cost) - uy * sint,
uy * uz * (1 - cost) + ux * sint,
(uz**2) * (1 - cost) + cost]])
return rotmat

def convert_for_lammps(atoms):
"""
Convert a parallel piped (forming right hand basis)
to lower triangular, low-tilt box that LAMMPS will accept.

This code draws from a previous LAMMPS interface:
https://svn.fysik.dtu.dk/projects/ase-extra/trunk/
ase-extra/trunk/ase/calculators/lammpslib.py
"""
ase_cell = atoms.get_cell()
cell = np.matrix.transpose(ase_cell)
# rotate bases into triangular matrix
tri_mat = np.zeros((3, 3))
A = cell[:, 0]
B = cell[:, 1]
C = cell[:, 2]
tri_mat[0, 0] = np.linalg.norm(A)
Ahat = A / np.linalg.norm(A)
AxBhat = np.cross(A, B) / np.linalg.norm(np.cross(A, B))
tri_mat[0, 1] = np.dot(B, Ahat)
tri_mat[1, 1] = np.linalg.norm(np.cross(Ahat, B))
tri_mat[0, 2] = np.dot(C, Ahat)
tri_mat[1, 2] = np.dot(C, np.cross(AxBhat, Ahat))
tri_mat[2, 2] = np.linalg.norm(np.dot(C, AxBhat))

atoms.set_cell(tri_mat.T, scale_atoms=True)
atoms.wrap(pbc=True)

# "flip" the cell if it is too skewed
newcell = atoms.get_cell()
while True:
xx, yy = newcell[0, 0], newcell[1, 1]
xy, xz, yz = newcell[1, 0], newcell[2, 0], newcell[2, 1]
cond1 = 2 * abs(xy) > xx
cond2 = 2 * abs(xz) > xx
cond3 = 2 * abs(yz) > yy
if not cond1 and not cond2 and not cond3:
break
if cond1:
newcell[1, 0] += xx * np.round((0.5 * xx - xy) / xx - 0.5)
if cond2:
newcell[2, 0] += xx * np.round((0.5 * xx - xz) / xx - 0.5)
if cond3:
newcell[2, 1] += yy * np.round((0.5 * yy - yz) / yy - 0.5)
newcell[2, 0] += xy * np.round((0.5 * yy - yz) / yy - 0.5)

atoms.set_cell(newcell, scale_atoms=False)
atoms.wrap(pbc=True)