"""
Classes and functions to create surface models.
Copyright Schrodinger, LLC. All rights reserved."""
import argparse
import numpy
from scipy import constants
from schrodinger import structure
from schrodinger.application.desmond.constants import IS_INFINITE
from schrodinger.application.matsci import clusterstruct
from schrodinger.application.matsci import msutils
from schrodinger.application.matsci import textlogger
from schrodinger.application.matsci.nano import particle
from schrodinger.application.matsci.nano import plane
from schrodinger.application.matsci.nano import xtal
from schrodinger.infra import mm
from schrodinger.infra import structure as infrastructure
from schrodinger.structutils import build
from schrodinger.structutils import color
from schrodinger.structutils import transform
XYZ_VECTORS = [
numpy.array(transform.X_AXIS),
numpy.array(transform.Y_AXIS),
numpy.array(transform.Z_AXIS)
]
OUTFILE_SUFFIX = '-surface'
H_INDEX_DEFAULT = 1
K_INDEX_DEFAULT = 1
L_INDEX_DEFAULT = 0
MIN_IDX_DEFAULT = 0
MAX_IDX_DEFAULT = 1
H_INDEX_MIN_DEFAULT = MIN_IDX_DEFAULT
H_INDEX_MAX_DEFAULT = MAX_IDX_DEFAULT
K_INDEX_MIN_DEFAULT = MIN_IDX_DEFAULT
K_INDEX_MAX_DEFAULT = MAX_IDX_DEFAULT
L_INDEX_MIN_DEFAULT = MIN_IDX_DEFAULT
L_INDEX_MAX_DEFAULT = MAX_IDX_DEFAULT
NORMAL_C_DEFAULT = False
BOTTOM_DEFAULT = 0.0
SLAB_THICKNESS_DEFAULT = 1.0
VACUUM_THICKNESS_DEFAULT = 2.0
H_INDEX_KEY = 'i_matsci_h_Miller_idx'
K_INDEX_KEY = 'i_matsci_k_Miller_idx'
L_INDEX_KEY = 'i_matsci_l_Miller_idx'
HKL_INDEX_KEYS = [H_INDEX_KEY, K_INDEX_KEY, L_INDEX_KEY]
TERMINAL_FRAGMENT_NONE = 'none'
TERMINAL_BOND_LENGTH = 1.2 # Ang.
OVERLAP_ATOM_THRESHOLD = 0.1 # Ang, different than xtal default (MATSCI-8902)
[docs]class ParserWrapper(object):
"""
Manages the argparse module to parse user command line
arguments.
"""
[docs] def __init__(self, scriptname, description):
"""
Create a ParserWrapper instance and process it.
:type scriptname: str
:param scriptname: name of this script
:type description: str
:param description: description of this script
"""
name = '$SCHRODINGER/run ' + scriptname
self.parserobj = argparse.ArgumentParser(
prog=name,
description=description,
add_help=False,
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
[docs] def loadIt(self):
"""
Load ParserWrapper with all options.
"""
self.loadHKLOptions()
self.loadOptions()
self.loadRequired()
self.loadCommon()
[docs] def loadHKLOptions(self):
"""
Load ParserWrapper with hkl options.
"""
h_index_help = """Specify the Miller integer h-index of the crystal
plane to be exposed in the surface model."""
self.parserobj.add_argument('-h_index',
type=int,
default=H_INDEX_DEFAULT,
help=h_index_help)
k_index_help = """Specify the Miller integer k-index of the crystal
plane to be exposed in the surface model."""
self.parserobj.add_argument('-k_index',
type=int,
default=K_INDEX_DEFAULT,
help=k_index_help)
l_index_help = """Specify the Miller integer l-index of the crystal
plane to be exposed in the surface model."""
self.parserobj.add_argument('-l_index',
type=int,
default=L_INDEX_DEFAULT,
help=l_index_help)
[docs] def loadOptions(self):
"""
Load ParserWrapper with options.
"""
self.parserobj.add_argument('-normal_c',
action='store_true',
help="Enforce normal C-axis.")
bottom_help = """Specify the bottom of the slab used for this surface
model. This option takes a distance along the plane normal in
fractional coordinates. Together with the -slab_thickness option an
envelope is defined for which any pair of planes may be exposed in
the surface model."""
self.parserobj.add_argument('-bottom',
type=float,
default=BOTTOM_DEFAULT,
help=bottom_help)
slab_thickness_help = """Specify the thickness of the slab used for this
surface model. This option takes a distance along the plane normal
in fractional coordinates. See also the -bottom option."""
self.parserobj.add_argument('-slab_thickness',
type=float,
default=SLAB_THICKNESS_DEFAULT,
help=slab_thickness_help)
vacuum_thickness_help = """Specify the thickness of the vacuum layer used
for this surface model. This option takes a distance along the plane
normal in fractional coordinates."""
self.parserobj.add_argument('-vacuum_thickness',
type=float,
default=VACUUM_THICKNESS_DEFAULT,
help=vacuum_thickness_help)
terminal_fragment_help = """Choose a fragment from the following list to
terminate the surface model: %s."""
self.parserobj.add_argument('-terminal_fragment',
type=str,
choices=particle.TERM_FRAGS,
default=TERMINAL_FRAGMENT_NONE,
help=terminal_fragment_help %
', '.join(particle.TERM_FRAGS))
[docs] def loadEnumeration(self):
"""
Load ParserWrapper with enumeration options.
"""
h_index_min_help = """Specify a starting Miller h-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-h_index_min',
type=int,
default=H_INDEX_MIN_DEFAULT,
help=h_index_min_help)
h_index_max_help = """Specify an ending Miller h-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-h_index_max',
type=int,
default=H_INDEX_MAX_DEFAULT,
help=h_index_max_help)
k_index_min_help = """Specify a starting Miller k-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-k_index_min',
type=int,
default=K_INDEX_MIN_DEFAULT,
help=k_index_min_help)
k_index_max_help = """Specify an ending Miller k-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-k_index_max',
type=int,
default=K_INDEX_MAX_DEFAULT,
help=k_index_max_help)
l_index_min_help = """Specify a starting Miller l-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-l_index_min',
type=int,
default=L_INDEX_MIN_DEFAULT,
help=l_index_min_help)
l_index_max_help = """Specify an ending Miller l-index to use
in enumerating surfaces."""
self.parserobj.add_argument('-l_index_max',
type=int,
default=L_INDEX_MAX_DEFAULT,
help=l_index_max_help)
hkl_indices_help = """Specify whitespace separated integers which
are internally grouped by order into triples of Miller indices
for sought surfaces."""
self.parserobj.add_argument('-hkl_indices',
type=int,
nargs='+',
help=hkl_indices_help)
[docs] def loadRequired(self):
"""
Load ParserWrapper with required options.
"""
input_file_help = """Specify a Maestro file containing a crystalline
cell from which to create the surface model."""
self.parserobj.add_argument('-input_file',
type=str,
required=True,
help=input_file_help)
[docs] def loadCommon(self):
"""
Load ParserWrapper with common options.
"""
self.parserobj.add_argument('-verbose',
action='store_true',
help='Turn on verbose printing.')
self.parserobj.add_argument('-help',
'-h',
action='help',
default=argparse.SUPPRESS,
help='Show this help message and exit.')
[docs] def parseArgs(self, args):
"""
Parse the command line arguments.
:type args: tuple
:param args: command line arguments
"""
self.options = self.parserobj.parse_args(args)
[docs]class Surface(object):
"""
Manage the building of a surface model.
"""
NUMDIGITS = 3
MSGWIDTH = 100
[docs] def __init__(self,
cell,
h_index=H_INDEX_DEFAULT,
k_index=K_INDEX_DEFAULT,
l_index=L_INDEX_DEFAULT,
normal_c=NORMAL_C_DEFAULT,
bottom=BOTTOM_DEFAULT,
slab_thickness=SLAB_THICKNESS_DEFAULT,
vacuum_thickness=VACUUM_THICKNESS_DEFAULT,
terminal_fragment=TERMINAL_FRAGMENT_NONE,
overlap_threshold=OVERLAP_ATOM_THRESHOLD,
do_bonding=None,
do_bond_orders=True,
logger=None):
"""
Create an instance.
:type cell: `schrodinger.structure.Structure`
:param cell: the crystalline cell used to define the surface model
:type h_index: int
:param h_index: the Miller h-index of the crystal plane to expose
:type k_index: int
:param k_index: the Miller k-index of the crystal plane to expose
:type l_index: int
:param l_index: the Miller l-index of the crystal plane to expose
:type bottom: float or None
:param bottom: the distance along the plane normal in fractional
coordinates to serve as the bottom of the surface model
:type slab_thickness: float or None
:param slab_thickness: thickness of the crystalline slab along the
plane normal in fractional coordinates
:type vacuum_thickness: float or None
:param vacuum_thickness: thickness of the vacuum layer along the
plane normal in fractional coordinates
:type terminal_fragment: str
:param terminal_fragment: the fragment by which to terminate the surface
:param float overlap_threshold: (Ang) distance used to define
overlapping atoms
:type do_bonding: bool or None
:param do_bonding: whether to do the bonding, i.e. connecting atoms,
see do_bond_orders kwarg for assigning bond orders, if None it is
determined from the other input
:type do_bond_orders: bool
:param do_bond_orders: whether to do the bond orders, used only if
bonding is enabled
:type logger: logging.Logger
:param logger: output logger
"""
self.cell = cell.copy()
self.normal_c = normal_c
self.bottom = bottom
self.slab_thickness = slab_thickness
self.vacuum_thickness = vacuum_thickness
assert terminal_fragment in particle.TERM_FRAGS
self.terminal_fragment = terminal_fragment
self.overlap_threshold = overlap_threshold
self.logger = logger
self.hkl = [h_index, k_index, l_index]
self.surface = None
self.is_infinite = xtal.is_infinite2(self.cell)
self.w_length = None
self._reduceHKL()
self._setDoBonding(do_bonding)
self.do_bond_orders = do_bond_orders
def _reduceHKL(self):
"""
Reduce HKLs to the smallest HKL indices and shift the cell accordingly.
Set reduced HKLs to self.reduced_hkl.
"""
# See also MATSCI-8146
self.reduced_hkl, divisor = plane.reduce_hkl(self.hkl)
# Nothing to reduce
if divisor == 1:
return
oshift = [0, 0, 0]
for idx, val in enumerate(self.reduced_hkl):
if val == 0:
continue
oshift[idx] = 1 / divisor
# Null shift (should never happen)
assert sum(oshift) != 0.
# Shift factional coordinates based on the divisor
vecs = xtal.get_vectors_from_chorus(self.cell)
fracs = xtal.trans_cart_to_frac_from_vecs(self.cell.getXYZ(), *vecs)
fracs_shift = fracs + oshift
self.cell.setXYZ(xtal.trans_frac_to_cart_from_vecs(fracs_shift, *vecs))
self.cell = xtal.move_atoms_into_cell(self.cell, preserve_bonding=True)
def _setDoBonding(self, do_bonding):
"""
Set the do bonding attribute.
:type do_bonding: bool or None
:param do_bonding: whether to do the bonding, i.e. both connecting
atoms and assigning bond orders, if None it is determined from
the other input
"""
assert self.is_infinite is not None
if self.terminal_fragment != TERMINAL_FRAGMENT_NONE:
self.do_bonding = True
return
if not self.is_infinite:
# Don't rebond molecular systems (see RB 52003)
self.do_bonding = False
return
if do_bonding is None:
space_group = self.cell.property.get(xtal.Crystal.SPACE_GROUP_KEY,
xtal.P1_SPACE_GROUP_SYMBOL)
if space_group != xtal.P1_SPACE_GROUP_SYMBOL:
self.do_bonding = True
return
ct_type = self.cell.property.get('s_ffio_ct_type')
if ct_type:
self.do_bonding = False
return
if (self.reduced_hkl in [[1, 0, 0], [0, 1, 0], [0, 0, 1]] and
self.slab_thickness == 1):
self.do_bonding = False
return
self.do_bonding = True
return
self.do_bonding = do_bonding
[docs] def getLatticeProperties(self, cell):
"""
Get the lattice properties.
:raise: ValueError if missing PDB PBC
"""
try:
lattice_parameters = xtal.get_lattice_param_properties(cell)
except KeyError:
keys = [
'A_KEY', 'B_KEY', 'C_KEY', 'ALPHA_KEY', 'BETA_KEY', 'GAMMA_KEY'
]
keys = [getattr(xtal.Crystal, key) for key in keys]
msg = ('The given cell is missing one or more of the '
'following properties: %s.') % keys
if self.logger:
self.logger.error(msg)
raise ValueError(msg)
a_vec, b_vec, c_vec = xtal.get_lattice_vectors(*lattice_parameters)
return list(lattice_parameters), [a_vec, b_vec, c_vec]
def _getBulkExtents(self):
"""
Return the bulk extents necessary for the slab build.
:rtype: list
:return: the bulk extents
"""
extents = [1, 1, int(numpy.ceil(self.bottom + self.slab_thickness))]
if self.is_infinite and self.terminal_fragment != TERMINAL_FRAGMENT_NONE:
extents[2] += 1
return extents
def _fixFiniteSystem(self):
"""
Given a surface with all atoms in the first cell, re-create the
molecules anchoring to the atoms with the smallest Z coordinate. After,
if any atoms are outside the cell, move the surface along Z by this
amount.
"""
anchors = []
xyz = self.surface.getXYZ()
for mol in self.surface.molecule:
min_z, anchor = numpy.Inf, None
for atom in mol.atom:
if atom.z < min_z:
min_z = atom.z
anchor = atom.index
anchors.append(anchor)
clusterstruct.contract_by_molecule(self.surface, anchors=anchors)
# Move the surface if needed
xyz = self.surface.getXYZ()
min_z = xyz[:, 2].min(axis=0)
if min_z < 0.:
transform.translate_structure(self.surface, 0., 0., -min_z)
[docs] def buildSlab(self):
"""
Build the slab.
"""
# Ensure that determinant is positive
if numpy.linalg.det(self.plane.basis) < 0.0:
if self.logger:
self.logger.info(
'Determinant is negative for the '
'transformation, will multiply the basis by -1.')
self.plane.basis *= -1
# Surface is in general not normal in the c-axis
self.surface = xtal.transform_pbc(self.cell,
self.plane.basis,
scale_only=False,
cell_only=False)
if self.do_bonding:
if self.do_bond_orders:
bond_orders = xtal.ParserWrapper.CHOICE_ON
else:
bond_orders = xtal.ParserWrapper.CHOICE_OFF
xtal_kwargs = {
'bonding': xtal.ParserWrapper.CHOICE_ON,
'bond_orders': bond_orders
}
else:
xtal_kwargs = {
'bonding': xtal.ParserWrapper.CHOICE_OFF,
'bond_orders': xtal.ParserWrapper.CHOICE_OFF
}
extents = self._getBulkExtents()
if self.is_infinite:
xtal_kwargs['translate'] = xtal.ParserWrapper.CHOICE_ON
else:
xtal_kwargs['translate'] = xtal.ParserWrapper.CHOICE_OFF
# Move entire molecules into the first cell (MATSCI-8357)
# See MATSCI-6963 and MATSCI-7935 where we want to ensure complete
# layers prior to adding vacuum
xtal_kwargs['translate_centroids'] = True
self.surface = xtal.get_cell(self.surface,
extents=extents,
xtal_kwargs=xtal_kwargs)
# Sync lattice params with chorus
xtal.set_pbc_properties(self.surface,
xtal.get_chorus_properties(self.surface))
# This `if` block was disabled in favor of
# xtal_kwargs['translate_centroids'] = True
# call above. See RB 52003 for the discussion
#if not self.is_infinite:
# xtal.translate_molecules(self.surface, fract_offset=0.1, maxes=True)
# handle normal cell, in order to provide a physically meaningful
# cell PBC bonds should be recalculated because the relative
# positioning of the top and bottom layers may have changed after angular
# strain, the get_normal_surf function only applies angular strain,
# not strain along the c-vector, and restores the original vacuum
# space at the top of the cell if it is infinite
if self.normal_c:
self.surface = xtal.get_normal_surf(
self.surface,
restore_vacuum=self.is_infinite,
overlap_threshold=self.overlap_threshold)
if self.do_bonding:
xtal.connect_atoms(self.surface)
if self.do_bond_orders:
self.surface = xtal.assign_bond_orders(self.surface)
if not self.is_infinite:
self._fixFiniteSystem()
self.surface = xtal.make_p1(self.surface)
vecs = numpy.array(xtal.get_vectors_from_chorus(self.surface))
self.w_length = transform.get_vector_magnitude(vecs[2] / extents[2])
def _getTerminalPairsDict(self, outside_idxs, inside_idxs, reverse=False):
"""
Return a dictionary containing bonding atom index pairs for
bonds to be terminated.
:type outside_idxs: list
:param outside_idxs: atom indices that are outside the cell,
sorted by increasing f3 coordinate
:type inside_idxs: set
:param inside_idxs: atom indices that are inside the cell
:type reverse: bool
:param reverse: if True traverse the outside atom indices in
reverse
:rtype: dict
:return: keys are outside atom indices, values are
sets of inside atom indices that are bonded to the given
outside atom
"""
# since below only differences in fractionals are considered
# the original vecs may be used, as opposed to the updated vecs
# from .pruneIt
vecs = xtal.get_vectors_from_chorus(self.surface)
outside_idxs = list(outside_idxs)
if reverse:
outside_idxs.reverse()
terminal_pairs_dict = {}
prev_f3 = None
for outside_idx in outside_idxs:
t_inside_idxs = set()
for neighbor in self.surface.atom[outside_idx].bonded_atoms:
if neighbor.index in inside_idxs:
t_inside_idxs.add(neighbor.index)
if t_inside_idxs:
terminal_pairs_dict[outside_idx] = t_inside_idxs
continue
# at this point no bonds to inside atoms were found for the
# current outside atom, because the traversal of outside atoms
# is sorted by f3 this would ordinarily mean to break out of the
# loop because you would be guaranteed that no other outside atoms
# would be bound to inside atoms, however before breaking out
# the entire f3 layer of atoms needs to be considered
curr_f3 = xtal.trans_cart_to_frac_from_vecs(
self.surface.atom[outside_idx].xyz, *vecs)[2]
if prev_f3 is None or abs(curr_f3 - prev_f3) <= xtal.FRACT_OFFSET:
prev_f3 = curr_f3
continue
break
return terminal_pairs_dict
def _makeOneToOneBonds(self, terminal_pairs_dict):
"""
Make one-to-many bonds one-to-one bonds.
:type terminal_pairs_dict: dict
:param terminal_pairs_dict: contains bonding pairs that need termination,
keys are outside atom indices, values are sets of inside atom indices
:rtype: dict
:return: contains bonding pairs that need termination, keys are outside
atom indices, values are inside atom indices
"""
# some outside indices are bonded to multiple inside indices so
# create unique bonds
maximally_bonded_idxs = set()
_terminal_pairs_dict = {}
for outside_idx, inside_idxs in terminal_pairs_dict.items():
outside_atom = self.surface.atom[outside_idx]
for inside_idx in inside_idxs:
inside_atom = self.surface.atom[inside_idx]
if inside_atom.bond_total == mm.MMCT_MAXBOND:
maximally_bonded_idxs.add(inside_idx)
continue
atom = self.surface.addAtom('H', *outside_atom.xyz)
self.surface.addBond(inside_idx, atom.index,
structure.BondType.Single)
_terminal_pairs_dict[atom.index] = inside_idx
amap = self.surface.deleteAtoms(terminal_pairs_dict.keys(),
renumber_map=True)
if maximally_bonded_idxs and self.logger:
maximally_bonded_idxs = [amap[idx] for idx in maximally_bonded_idxs]
msg = ('The following slab atoms are maximally bonded with '
f'{mm.MMCT_MAXBOND} bonds: {maximally_bonded_idxs}.')
self.logger.warning(msg)
_terminal_pairs_dict = {
amap[k]: amap[v] for k, v in _terminal_pairs_dict.items()
}
return _terminal_pairs_dict
def _makeConventionalBonds(self, terminal_pairs_dict):
"""
Make conventional bonds from bonds that cross the PBC.
:type terminal_pairs_dict: dict
:param terminal_pairs_dict: contains bonding pairs that need termination,
keys are outside atom indices, values are inside atom indices
"""
pbc = infrastructure.PBC(self.surface)
for outside_idx, inside_idx in terminal_pairs_dict.items():
outside_atom = self.surface.atom[outside_idx]
inside_atom = self.surface.atom[inside_idx]
outside_atom.xyz = pbc.getNearestImage(self.surface, inside_atom,
outside_atom)
self.surface.adjust(TERMINAL_BOND_LENGTH, inside_idx, outside_idx)
[docs] def pruneIt(self):
"""
Prune the slab cell.
:raise: ValueError if the slab model has zero atoms
:rtype: dict
:return: contains bonding pairs that need termination,
keys are outside atom indices, values are inside atom indices
"""
vecs = numpy.array(xtal.get_vectors_from_chorus(self.surface))
fracs = xtal.trans_cart_to_frac_from_vecs(self.surface.getXYZ(), *vecs)
# bottom needs to be updated to correspond with the extended cell
extents = self._getBulkExtents()
bottom = self.bottom / extents[2]
# shift cartesians to define the slab bottom
fracs[:, 2] -= bottom
self.surface.setXYZ(xtal.trans_frac_to_cart_from_vecs(fracs, *vecs))
# define final w-vector (or c-vector)
vecs[2] = transform.get_normalized_vector(vecs[2]) * self.slab_thickness * \
self.w_length
# collect atom indices below and above the final cell
fracs = xtal.trans_cart_to_frac_from_vecs(self.surface.getXYZ(), *vecs)
outside_below_idxs = list(
numpy.where(fracs[:, 2] < -xtal.FRACT_OFFSET)[0] + 1)
outside_above_idxs = list(
numpy.where(fracs[:, 2] >= 1 - xtal.FRACT_OFFSET)[0] + 1)
outside_idxs = outside_below_idxs + outside_above_idxs
# molecules that cross the cell boundary should be complete
if not self.is_infinite:
outside_idxs = set(outside_idxs)
for mol in self.surface.molecule:
all_idxs = set(mol.getAtomIndices())
if not all_idxs.issubset(outside_idxs):
outside_idxs.difference_update(all_idxs)
if len(outside_idxs) == self.surface.atom_total:
msg = ('The specified slab bottom and thickness result in a model '
'with zero atoms.')
raise ValueError(msg)
if not self.is_infinite or self.terminal_fragment == TERMINAL_FRAGMENT_NONE:
xtal.set_pbc_properties(self.surface, vecs.flat)
self.surface.deleteAtoms(outside_idxs)
return {}
# build a dictionary containing all atom index pairs that need
# termination, keys are outside indices, values are sets of inside
# atom indices
outside_below_idxs = sorted(outside_below_idxs,
key=lambda x: fracs[x - 1, 2])
outside_above_idxs = sorted(outside_above_idxs,
key=lambda x: fracs[x - 1, 2])
inside_idxs = set(range(1, self.surface.atom_total +
1)).difference(outside_idxs)
if outside_below_idxs:
terminal_pairs_dict_below = self._getTerminalPairsDict(
outside_below_idxs, inside_idxs, reverse=True)
else:
terminal_pairs_dict_below = self._getTerminalPairsDict(
outside_above_idxs, inside_idxs, reverse=True)
terminal_pairs_dict_above = self._getTerminalPairsDict(
outside_above_idxs, inside_idxs, reverse=False)
terminal_pairs_dict = terminal_pairs_dict_below.copy()
for k, v in terminal_pairs_dict_above.items():
terminal_pairs_dict.setdefault(k, set()).update(v)
# remove outside indices that not not bonded to inside indices
to_delete = set(outside_idxs).difference(terminal_pairs_dict)
if to_delete:
amap = self.surface.deleteAtoms(to_delete, renumber_map=True)
terminal_pairs_dict = {
amap[k]: list(map(lambda x: amap[x], vs))
for k, vs in terminal_pairs_dict.items()
}
terminal_pairs_dict = self._makeOneToOneBonds(terminal_pairs_dict)
self._makeConventionalBonds(terminal_pairs_dict)
xtal.set_pbc_properties(self.surface, vecs.flat)
return terminal_pairs_dict
[docs] def doTermination(self, terminal_pairs_dict):
"""
Do the termination.
:type terminal_pairs_dict: dict
:param terminal_pairs_dict: contains bonding pairs that need termination,
keys are outside atom indices, values are inside atom indices
"""
# terminate if not H
if self.terminal_fragment != particle.TERM_FRAG:
frag_group, frag_name = particle.TERM_DICT[self.terminal_fragment]
outside_idxs = list(terminal_pairs_dict.keys())
inside_idxs = list(terminal_pairs_dict.values())
while outside_idxs and inside_idxs:
outside_idx = outside_idxs.pop(0)
inside_idx = inside_idxs.pop(0)
amap = build.attach_fragment(self.surface, inside_idx,
outside_idx, frag_group, frag_name,
particle.TERM_DIRECTION)
outside_idxs = list(map(lambda x: amap[x], outside_idxs))
inside_idxs = list(map(lambda x: amap[x], inside_idxs))
vecs = numpy.array(xtal.get_vectors_from_chorus(self.surface))
fracs = xtal.trans_cart_to_frac_from_vecs(self.surface.getXYZ(), *vecs)
# update the Cartesian coordinates of the terminated slab so that
# the bottom terminating atoms are inside the cell
fracs[:, 2] += abs(min(fracs[:, 2]))
self.surface.setXYZ(xtal.trans_frac_to_cart_from_vecs(fracs, *vecs))
# expand the top of the PBC to include the terminating
# atoms plus a vacuum buffer of 2 Ang.
vecs[2] *= max(fracs[:, 2])
vecs[2] += 2 * transform.get_normalized_vector(vecs[2])
xtal.set_pbc_properties(self.surface, vecs.flat)
color.apply_color_scheme(self.surface, 'element')
[docs] def addVacuum(self):
"""
Add vacuum.
"""
vecs = numpy.array(xtal.get_vectors_from_chorus(self.surface))
# for infinite non-terminated slabs with integer thickness remove bonds
# that cross the PBC along the c-vector now that vacuum is being added
if self.is_infinite and self.terminal_fragment == TERMINAL_FRAGMENT_NONE and \
not self.bottom and not numpy.fmod(self.slab_thickness, 1):
xtal.delete_pbc_bonds(self.surface, along_vector_idxs=[2])
vecs[2] = transform.get_normalized_vector(vecs[2]) * (
self.slab_thickness + self.vacuum_thickness) * self.w_length
xtal.set_pbc_properties(self.surface, vecs.flat)
[docs] def setSurfaceProperties(self):
"""
Set surface properties.
"""
xtal.set_physical_properties(self.surface)
xtal.clear_asu_and_fractional_properties(self.surface)
# atomic properties
crystal_keys = [xtal.Crystal.SYMMETRY_LABEL_KEY, xtal.COV_RADIUS_KEY]
for key in crystal_keys:
msutils.remove_atom_property(self.surface, key)
# set title
miller_indices = ''.join(map(str, self.hkl))
self.surface.title = self.cell.title + '-' + miller_indices + OUTFILE_SUFFIX
# set Miller idxs
for key, value in zip(HKL_INDEX_KEYS, self.hkl):
self.surface.property[key] = value
self.surface.property[IS_INFINITE] = self.is_infinite
self.surface.property[xtal.PBC_POSITION_KEY] = \
xtal.ANCHOR_PBC_POSITION % ('0', '0', '0')
[docs] def printParams(self):
"""
Log the parameters.
"""
lparams, lvecs = self.getLatticeProperties(self.cell)
sparams, svecs = self.getLatticeProperties(self.surface)
HEADER = 'Surface Parameters'
lattice_parameters = \
textlogger.get_param_string('Lattice (a, b, c, alpha, beta, gamma) / Ang. & degree',
', '.join([str(round(param, self.NUMDIGITS)) for param in lparams]),
self.MSGWIDTH)
slab_parameters = \
textlogger.get_param_string('Slab (a, b, c, alpha, beta, gamma) / Ang. & degree',
', '.join([str(round(param, self.NUMDIGITS)) for param in sparams]),
self.MSGWIDTH)
miller_indices = ''.join(map(str, self.hkl))
miller_indices = textlogger.get_param_string('Miller indices',
miller_indices,
self.MSGWIDTH)
reduced_miller_indices = ''.join(map(str, self.reduced_hkl))
reduced_miller_indices = textlogger.get_param_string(
'Reduced Miller indices', reduced_miller_indices, self.MSGWIDTH)
inter_planar_separation = \
textlogger.get_param_string('Inter-planar separation / Ang.',
round(self.plane.inter_planar_separation, self.NUMDIGITS), self.MSGWIDTH)
bottom = textlogger.get_param_string('Bottom / n',
round(self.bottom,
2), self.MSGWIDTH)
slab_thickness = textlogger.get_param_string(
'Slab thickness / n', round(self.slab_thickness, 2), self.MSGWIDTH)
vacuum_thickness = textlogger.get_param_string(
'Vacuum thickness / n', round(self.vacuum_thickness, 2),
self.MSGWIDTH)
terminal_fragment = textlogger.get_param_string('Terminal fragment',
self.terminal_fragment,
self.MSGWIDTH)
a_param = textlogger.get_param_string('Length slab vector a / Ang.', \
round(sparams[0], self.NUMDIGITS), self.MSGWIDTH)
b_param = textlogger.get_param_string('Length slab vector b / Ang.', \
round(sparams[1], self.NUMDIGITS), self.MSGWIDTH)
c_param = textlogger.get_param_string('Length slab vector c / Ang.', \
round(sparams[2], self.NUMDIGITS), self.MSGWIDTH)
alpha_param = textlogger.get_param_string('Angle slab vectors b-c / Degree', \
round(sparams[3], self.NUMDIGITS), self.MSGWIDTH)
beta_param = textlogger.get_param_string('Angle slab vectors a-c / Degree', \
round(sparams[4], self.NUMDIGITS), self.MSGWIDTH)
gamma_param = textlogger.get_param_string('Angle slab vectors a-b / Degree', \
round(sparams[5], self.NUMDIGITS), self.MSGWIDTH)
u_vec = textlogger.get_param_string('Slab vector a (planar) / Ang.', \
svecs[0], self.MSGWIDTH)
v_vec = textlogger.get_param_string('Slab vector b (planar) / Ang.', \
svecs[1], self.MSGWIDTH)
w_vec = textlogger.get_param_string('Slab vector c (thickness) / Ang.', \
svecs[2], self.MSGWIDTH)
self.logger.info(HEADER)
self.logger.info('-' * len(HEADER))
self.logger.info(lattice_parameters)
self.logger.info(slab_parameters)
self.logger.info(miller_indices)
self.logger.info(reduced_miller_indices)
self.logger.info(inter_planar_separation)
self.logger.info(bottom)
self.logger.info(slab_thickness)
self.logger.info(vacuum_thickness)
self.logger.info(terminal_fragment)
self.logger.info(a_param)
self.logger.info(b_param)
self.logger.info(c_param)
self.logger.info(alpha_param)
self.logger.info(beta_param)
self.logger.info(gamma_param)
self.logger.info(u_vec)
self.logger.info(v_vec)
self.logger.info(w_vec)
self.logger.info('')
[docs] def runIt(self):
"""
Create the surface model.
"""
lparams, lvecs = self.getLatticeProperties(self.cell)
self.plane = plane.CrystalPlane(*self.reduced_hkl, *lvecs)
self.plane.getSlabVectors()
self.buildSlab()
terminal_pairs_dict = self.pruneIt()
if terminal_pairs_dict:
self.doTermination(terminal_pairs_dict)
if self.vacuum_thickness > 0:
self.addVacuum()
self.setSurfaceProperties()
if self.logger:
self.printParams()
[docs]def get_hkl(astructure):
"""
Return an hkl Miller index triple for the given structure.
If any of the structure properties are missing then by
convention return (0, 0, 1).
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure for which you want the Miller
indices
:rtype: tuple
:return: hkl Miller index triple
"""
try:
hkl = tuple([astructure.property[key] for key in HKL_INDEX_KEYS])
except KeyError:
hkl = (0, 0, 1)
return hkl
[docs]def set_hkl(astructure, hkl):
"""
Set hkl Miller indices to a given structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure for which you want the Miller
indices
:type hkl: list
:param hkl: hkl Miller indices
"""
assert len(hkl) == len(HKL_INDEX_KEYS)
for prop, idx in zip(HKL_INDEX_KEYS, hkl):
astructure.property[prop] = idx
[docs]def maestro_rotate_cell(astructure, origin, a_vec, b_vec, c_vec):
"""
Rotate the cell for proper Maestro view.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure that you want rotated
:type origin: numpy.array
:param origin: the origin
:type a_vec: numpy.array
:param a_vec: the a lattice vector
:type b_vec: numpy.array
:param b_vec: the b lattice vector
:type c_vec: numpy.array
:param c_vec: the c lattice vector
:rtype: `schrodinger.structure.Structure` and four numpy.array
:return: rotated structure and vectors
"""
# the structure and lattice vectors need to be rotated such that
# the a-vec is aligned along the x-axis and the c-vec is aligned
# along the z-axis unless the gamma angle, i.e. the angle between
# the a-vec and b-vec, is right in which case the b-vec should
# be aligned with the y-axis
def rotate(index, vectors):
rotation = transform.get_alignment_matrix(vectors[index],
XYZ_VECTORS[index])
rotated_vectors = []
for vector in vectors:
rotated_vectors.append(
transform.transform_atom_coordinates(vector, rotation))
transform.transform_structure(astructure, rotation)
return rotated_vectors
trans_vec = numpy.array(xtal.ParserWrapper.ORIGIN) - origin
transform.translate_structure(astructure, *trans_vec)
gamma = transform.get_angle_between_vectors(a_vec, b_vec) / constants.degree
tmp_vectors = rotate(0, [a_vec, b_vec, c_vec])
index = 1
if gamma != 90.0:
index = 2
a_vec_p, b_vec_p, c_vec_p = rotate(index, tmp_vectors)
return astructure, numpy.array(
xtal.ParserWrapper.ORIGIN), a_vec_p, b_vec_p, c_vec_p