"""
Classes and functions to create nanoparticles.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
import argparse
import itertools
import math
import warnings
from collections import OrderedDict
from past.utils import old_div
import numpy
from scipy import constants as scipy_constants
from schrodinger.application.matsci import shapes
from schrodinger.application.matsci import textlogger
from schrodinger.application.matsci.nano import constants
from schrodinger.application.matsci.nano import plane
from schrodinger.application.matsci.nano import xtal
from schrodinger.structutils import build
from schrodinger.structutils import color
from schrodinger.structutils import transform
_version = '$Revision 0.0 $'
DEFAULT_SHAPE = shapes.Cube.NAME
ALL_SHAPES = list(shapes.SHAPES_NAMES_TO_OBJECTS_DICT)
OUTFILE_SUFFIX = '-out'
FRAGMENT_NAME = constants.Constants.FRAGMENT_NAME
TERM_FRAG = constants.Constants.TERMFRAG
TERM_FRAGS = constants.Constants.TERMFRAGS
TERM_DICT = constants.Constants.TERMDICT
TERM_DIRECTION = 'forward'
TERM_BOND_LENGTH = 1.50
NANOPARTICLE_NATOMS_KEY = 'i_matsci_Nanoparticle_Natoms'
NANOPARTICLE_VOLUME_KEY = 'r_matsci_Volume/Ang.^3'
NANOPARTICLE_SURFACE_AREA_KEY = 'r_matsci_Surface_Area/Ang.^2'
NANOPARTICLE_CIRCUMRADIUS_KEY = 'r_matsci_Circumradius/Ang.'
NANOPARTICLE_DENSITY_KEY = 'r_matsci_Density/g/cm^3'
NANOPARTICLE_TEMPLATE = 'b_matsci_Nanoparticle_Template'
ORIGINAL_ATOM_KEY = 'b_matsci_np_gui_original_atom'
NUM_DECIMAL = 4
CM_TO_ANGSTROM = float(100000000)
TOKEN_SEPARATOR = ','
KEY_VALUE_SEPARATOR = '='
AXIS_KEY = 'axis'
VECTOR_KEY = 'vector'
BASIS_KEY = 'basis'
HKL_VALUE = 'hkl'
ABC_VALUE = 'abc'
XYZ_VALUE = 'xyz'
ALIGNMENT_KEYS = [AXIS_KEY, VECTOR_KEY, BASIS_KEY]
BASIS_VALUES = [HKL_VALUE, ABC_VALUE, XYZ_VALUE]
PRIMARY_ALIGNMENT_DEFAULT = OrderedDict([(AXIS_KEY, '1'), (VECTOR_KEY, '1.0 0.0 0.0'), \
(BASIS_KEY, HKL_VALUE)])
SECONDARY_ALIGNMENT_DEFAULT = OrderedDict([(AXIS_KEY, '1'), (VECTOR_KEY, '0.0 1.0 0.0'), \
(BASIS_KEY, ABC_VALUE)])
dict_to_str = lambda x: [
TOKEN_SEPARATOR.join(
[KEY_VALUE_SEPARATOR.join(map(str, pair)) for pair in x.items()])
]
PRIMARY_ALIGNMENT_DEFAULT_LIST = dict_to_str(PRIMARY_ALIGNMENT_DEFAULT)
SECONDARY_ALIGNMENT_DEFAULT_LIST = dict_to_str(SECONDARY_ALIGNMENT_DEFAULT)
INCLUDE_AXES_DEFAULT = False
VERTICES_DEFAULT = []
ALLOW_FRAGMENTS = False
[docs]def get_lattice_properties():
"""
Get the properties which are required for defining a crystal structure
:rtype: list of str
:return: Each item in the list is a structure property that is required to
be present if a crystal structure is defined
"""
aclass = xtal.Crystal
return [
aclass.A_KEY, aclass.B_KEY, aclass.C_KEY, aclass.ALPHA_KEY,
aclass.BETA_KEY, aclass.GAMMA_KEY
]
[docs]def check_if_lattice_properties_exist(struct, keys=None):
"""
Check if the structure has all the required crystal lattice properties
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:type keys: list
:param keys: The list of crystal properties to check. If not given, the
default set of crystal properties from get_lattice_properties will be
checked
:rtype: bool
:return: If all the requested properties are present on the structure
"""
if not keys:
keys = get_lattice_properties()
return set(keys) <= set(struct.property)
[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 options.
"""
shape_help = """Specify the shape of the nanoparticle."""
self.parserobj.add_argument('-shape',
type=str,
choices=ALL_SHAPES,
default=DEFAULT_SHAPE,
help=shape_help)
origin_help = """When cutting the specified shape of nanoparticle from
the crystalline solid it is necessary to specify the placement of
the shape in the solid. This origin specifies the location of the
shape's centroid. See the options -primary_alignment and
-secondary_alignment which are also useful for orienting the shape.
Note that if the shape is being aligned along a crystal plane
(see -primary_alignment) the specified origin will be updated by
translating it along the plane normal."""
self.parserobj.add_argument('-origin',
type=float,
nargs=3,
default=shapes.ORIGIN,
metavar='float',
help=origin_help)
params_help = """Specify the set of parameters that define the geometry
of the shape. Units of Angstrom and degree are used. Here is a
collection of shapes, their default parameters, and parameter
descriptions: %s."""
defaults = []
for name in ALL_SHAPES:
values = []
shape_obj = shapes.get_shape_object_by_name(name)
for value, info in zip(shape_obj.DEFAULT, shape_obj.DESCRIPTION):
values.append(str(value) + ' (' + info + ')')
defaults.append(name + ': ' + ' '.join(values))
defaults = ', '.join(defaults)
self.parserobj.add_argument('-params',
type=float,
nargs='*',
metavar='float',
help=params_help % defaults)
vertices_help = """Overwrite the vertices of the polyhedron that has been
specified via the -origin, -params, -primary_alignment, and
-secondary_alignment options to this script. This offers the user
the most flexible level of control with regard to orienting the polyhedron
relative to the supercell. The polyhedron returned from this script
when run without this command can be imported into Maestro and
oriented relative to the supercell. The Cartesian coordinates
of the vertices can then be fed in here using the form: x_1 y_1
z_1 x_2 ... x_N y_N z_N and will be parsed three numbers at a
time. When using this option it is important that none of the
internal degrees of freedom that define the polyhedron be altered.
Avoid using this option with output from a run that uses the
-include_axes option."""
self.parserobj.add_argument('-vertices',
type=float,
nargs='*',
default=VERTICES_DEFAULT,
help=vertices_help)
primary_alignment_help = """Specify the primary alignment using the
following sub-options, %s. Sub-options should be separated from their
values using an '%s', while individual sub-option pairs should be
separated from each other using a '%s'. The %s sub-option takes an
integer axis ID that specifies which of the available symmetry unique
primary axes to use in the alignment. Primary axes are shape face
normals and are sorted by decreasing face area. Here is a collection
of shapes and the number of available primary axes: %s. The %s sub-option
takes three whitespace separated floating point numbers defining a
vector in a basis given by the %s sub-option. Available bases include
%s, which is reciprocal lattice space or the normal to the crystal
plane given by Miller indices h, k, and l, %s, which is direct lattice
space, and %s, which is Cartesian space. With this option the shape
will be rotated so that the normal of the face specified via the
primary axis ID lies parallel to the specified vector. If a crystal
plane has been specified then the face will be translated so that
it is coplanar with the crystal plane. As a consequence of this the
origin, specified with the origin option, will be updated. By default
the face of largest area, i.e. %s, will be aligned along the %s
crystal plane. Together with the -secondary alignment options the rigid
body rotations of the shape are fixed."""
values = format_alignment_ids('NUM_UNIQUE_FACES')
self.parserobj.add_argument(
'-primary_alignment',
type=str,
nargs='+',
default=PRIMARY_ALIGNMENT_DEFAULT_LIST,
help=primary_alignment_help %
(', '.join(ALIGNMENT_KEYS), KEY_VALUE_SEPARATOR, TOKEN_SEPARATOR,
AXIS_KEY, values, VECTOR_KEY, BASIS_KEY, HKL_VALUE, ABC_VALUE,
XYZ_VALUE, PRIMARY_ALIGNMENT_DEFAULT[AXIS_KEY],
PRIMARY_ALIGNMENT_DEFAULT[VECTOR_KEY]))
secondary_alignment_help = """Specify the secondary alignment using the
following sub-options, %s. Sub-options should be separated from their
values using an '%s', while individual sub-option pairs should be
separated from each other using a '%s'. The %s sub-option takes an
integer axis ID that specifies which of the available symmetry unique
secondary axes to use in the alignment. Secondary axes are shape face
edges of the primary face, specified with -primary_alignment, and are
sorted by decreasing edge length. Here is a collection of shapes and the
number of available secondary axes: %s. The %s sub-option takes three
whitespace separated floating point numbers defining a vector in a
basis given by the %s sub-option. Available bases include %s, which
is reciprocal lattice space or the normal to the crystal plane given
by Miller indices h, k, and l, %s, which is direct lattice space, and
%s, which is Cartesian space. With this option the shape will be
rotated, following the primary alignment, so that the edge, specified
via the secondary axis ID, of the primary face lies parallel to the
specified vector. If a crystal plane has been specified then the edge
will lie perpendicular to that plane. For convenience if the vectors
specified in the primary and secondary alignment are non-orthogonal
then rather than using this vector directly we will use it's component
that is orthogonal to the vector used in the primary alignment. By
default the longest edge, i.e. %s, of the face with the largest area
will be aligned along the %s direct lattice vector. Together with
the -primary alignment options the rigid body rotations of the shape
are fixed."""
values = format_alignment_ids('NUM_UNIQUE_EDGES')
self.parserobj.add_argument(
'-secondary_alignment',
type=str,
nargs='+',
default=SECONDARY_ALIGNMENT_DEFAULT_LIST,
help=secondary_alignment_help %
(', '.join(ALIGNMENT_KEYS), KEY_VALUE_SEPARATOR, TOKEN_SEPARATOR,
AXIS_KEY, values, VECTOR_KEY, BASIS_KEY, HKL_VALUE, ABC_VALUE,
XYZ_VALUE, SECONDARY_ALIGNMENT_DEFAULT[AXIS_KEY],
SECONDARY_ALIGNMENT_DEFAULT[VECTOR_KEY]))
include_axes_help = """Specifies that unit vectors defining the primary and
seconary alignment axes be added to the returned template structure of
the shape."""
self.parserobj.add_argument('-include_axes',
action='store_true',
help=include_axes_help)
term_frag_help = """Choose a fragment from the following list to
terminate the nanoparticle: %s. See the -allow_fragments option
for additional information."""
self.parserobj.add_argument('-term_frag',
type=str,
choices=TERM_FRAGS,
metavar=FRAGMENT_NAME,
default=TERM_FRAG,
help=term_frag_help % ', '.join(TERM_FRAGS))
term_bond_length_help = """Specify the bond length (in Angstrom) to use
when terminating the nanoparticle with fragments. This is the length
of the bond connecting the nanoparticle and fragment."""
self.parserobj.add_argument('-term_bond_length',
type=float,
default=TERM_BOND_LENGTH,
help=term_bond_length_help)
allow_fragments_help = """By default if the nanoparticle is being cut
from a molecular solid then molecules that lie across the shape
boundary are included in the nanoparticle. In this case the
-term_frag option is ignored and no termination is performed.
Use this option to instead fragment the molecule at the boundary
of the shape. In this case the -term_frag option is honored.
Note that because sometimes non-molecular solids may be incorrectly
classified as molecular solids, this option may sometimes be
necessary to force the cutting of nanoparticles from non-molecular
solids."""
self.parserobj.add_argument('-allow_fragments',
action='store_true',
help=allow_fragments_help)
input_file_help = """Specify a Maestro file containing a crystalline
supercell from which to cut the specified shape of nanoparticle."""
self.parserobj.add_argument('-input_file',
type=str,
required=True,
help=input_file_help)
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.')
self.parserobj.add_argument(
'-version',
'-v',
action='version',
default=argparse.SUPPRESS,
version=_version,
help="""Show the script's version number 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]def group_subtuples_by_key_and_value(atuple):
"""
Take a tuple of two-element subtuples and group them by both key
and value and collect into subsubtuples. For example, take
((1, 2), (1, 3), (4, 5), (6, 7), (8, 7)) and return (((1), (2, 3)),
((4), (5)), ((6, 8), (7))).
:type atuple: tuple
:param atuple: contains two-element tuples
:rtype: tuple
:return: contains two-element tuples each of which contains a keys
tuple and a values tuple
"""
atuple_grouped = ()
# sort by value
atuple_sorted = tuple(sorted(atuple, key=lambda x: x[1]))
# handle by value
next_atuple = []
for key, group in itertools.groupby(atuple_sorted, lambda x: x[1]):
things = sorted([thing[0] for thing in group])
if len(things) > 1:
atuple_grouped += ((tuple(things), (key,)),)
else:
next_atuple.append((things[0], key))
# sort by key
next_atuple_sorted = tuple(sorted(next_atuple, key=lambda x: x[0]))
# handle by key
for key, group in itertools.groupby(next_atuple_sorted, lambda x: x[0]):
things = sorted([thing[1] for thing in group])
atuple_grouped += (((key,), tuple(things)),)
return atuple_grouped
[docs]class Nanoparticle(object):
"""
Manage the building of a nanoparticle.
"""
MSGWIDTH = 100
[docs] def __init__(self,
supercell,
shape=DEFAULT_SHAPE,
origin=None,
params=None,
overwrite_vertices=VERTICES_DEFAULT,
primary_alignment_dict=PRIMARY_ALIGNMENT_DEFAULT,
secondary_alignment_dict=SECONDARY_ALIGNMENT_DEFAULT,
include_axes=INCLUDE_AXES_DEFAULT,
term_frag=TERM_FRAG,
term_bond_length=TERM_BOND_LENGTH,
allow_fragments=ALLOW_FRAGMENTS,
logger=None,
no_term_faces=None):
"""
Create an instance.
:type supercell: `schrodinger.structure.Structure`
:param supercell: the crystalline supercell from which to cut the
nanaparticle
:type shape: str
:param shape: the shape of the nanoparticle
:type origin: list
:param origin: list of three floats specifying the origin of
the shape in the crystal
:type params: list
:param params: list of parameters defining the specified shape
:type overwrite_vertices: list
:param overwrite_vertices: vertices to overwrite those of the created
shape, just of list of floats
:type primary_alignment_dict: dict
:param primary_alignment_dict: a dictionary defining the primary alignment
:type secondary_alignment_dict: dict
:param secondary_alignment_dict: a dictionary defining the secondary alignment
:type include_axes: bool
:param include_axes: include unit vectors for the primary and secondary alignment
axes in the template
:type term_frag: str
:param term_frag: the fragment by which to terminate the nanoparticle
:type term_bond_length: float
:param term_bond_length: the bond length of the bond connecting the
nanoparticle with the fragment
:type allow_fragments: bool
:param allow_fragments: specify if molecular fragments are allowed at the
boundary of the nanoparticle
:type logger: logging.Logger
:param logger: output logger
:type no_term_faces: list
:param no_term_faces: by default when terminating polyhedral nanoparticles all
faces are considered, this option specifies that the given integer faces not
be terminated (integers start at zero)
"""
if origin is None:
origin = shapes.ORIGIN
if not params:
params = shapes.get_shape_object_by_name(shape).DEFAULT
self.supercell = supercell
self.shape = shape
self.origin = numpy.array(origin)
self.params = params
self.overwrite_vertices = overwrite_vertices
self.primary_alignment_dict = primary_alignment_dict
self.secondary_alignment_dict = secondary_alignment_dict
self.include_axes = include_axes
self.term_frag = term_frag
self.term_bond_length = term_bond_length
self.allow_fragments = allow_fragments
self.logger = logger
# the face indicies must correspond with those given in
# shapes.<ShapeObject>.INDICES
if no_term_faces is None:
self.no_term_faces = []
else:
self.no_term_faces = no_term_faces
self.is_infinite = None
self.molecule_completion = None
self.lattice_parameters = None
self.plane = None
self.shape_object = None
self.template = None
self.inside_indices = None
self.boundary_indices = None
self.outside_indices = None
self.no_term_indices = set()
self.nanoparticle = self.supercell.copy()
self.polyhedron = None
self.nanoparticle_natoms = None
self.nanoparticle_volume = None
self.nanoparticle_surface_area = None
self.nanoparticle_circumradius = None
self.nanoparticle_density = None
self.embedded = None
[docs] def checkInput(self):
"""
Check user input.
"""
CheckInput().checkShapeParams(self.shape, self.params, self.logger)
self.primary_alignment_dict, self.secondary_alignment_dict = \
CheckInput().checkAlignment(self.shape, self.primary_alignment_dict,
self.secondary_alignment_dict, self.logger)
[docs] def setMoleculeCompletion(self):
"""
Set whether or not molecules should be completed at the
nanoparticle boundary.
"""
self.molecule_completion = not self.allow_fragments
if self.molecule_completion:
try:
self.is_infinite = xtal.is_infinite(self.nanoparticle)
except ValueError:
self.is_infinite = False
self.molecule_completion = not self.is_infinite
if not self.nanoparticle.bond:
self.molecule_completion = True
if self.molecule_completion:
self.term_frag = TERM_FRAGS[0]
[docs] def isInside(self, aatom):
"""
Return whether or not the given atom is inside the shape. If the
given atom has not yet been categorized then categorize it and
append the result to the collections.
:type aatom: `schrodinger.structure._StructureAtom`
:param aatom: the atom in question
:rtype: bool
:return: True if the atom is inside the shape, False otherwise
"""
if aatom.index in self.inside_indices:
inside = True
elif aatom.index in self.outside_indices:
inside = False
else:
inside = self.shape_object.pointInside(numpy.array(aatom.xyz))
if inside:
self.inside_indices.add(aatom.index)
else:
self.outside_indices.add(aatom.index)
# if we will be performing a termination then store
# this outside atom's index unless the
# atom-to-polyhedron-center line segment intersects
# a plane of a polyhedron face (not the face itself
# but rather the plane that the face resides in) that
# is supposed to be terminated
if self.polyhedron and self.term_frag != TERM_FRAGS[0] and \
self.no_term_faces:
terminatable = False
for aface in self.shape_object.faces:
if aface.index - 1 not in self.no_term_faces and \
aface.plane_intersected:
terminatable = True
break
if not terminatable:
self.no_term_indices.add(aatom.index)
return inside
[docs] def checkIntersectingBond(self, atom1, atom2):
"""
Check if the bond given by the specified atoms
can be intersected by the shape of nanoparticle.
We will only be either keeping a single atom from
this bond (no termination) or using it to terminate
with a fragment that is expected to connect via
a single bond.
:type atom1: int
:param atom1: the first atom index of the bond
:type atom2: int
:param atom2: the second atom index of the bond
"""
CAN_NOT_INTERSECT_MSG = (
'The nanoparticle specifications '
'that you have provided result in a shape that is '
'intersecting bonds in the crystalline supercell '
'that ought not to be intersected, i.e. those '
'of bond order greater than one. The nanoparticle '
'should be cut from the supercell at single bonds.')
if self.nanoparticle.getBond(atom1, atom2).order > 1:
if self.logger:
self.logger.error(CAN_NOT_INTERSECT_MSG)
raise ValueError(CAN_NOT_INTERSECT_MSG)
[docs] def partitionAtoms(self):
"""
Return three non-overlapping lists of atom indices, (1) for
atoms inside the shape, (2) for atoms belonging to bonds which
intersect the shape boundary (a (inside, outside) tuple of
atom indices for the intersecting bonds), and (3) for atoms
that lie outside of shape.
"""
self.inside_indices, self.boundary_indices, self.outside_indices = \
set(), set(), set()
# handle two cases of nanoparticles (1) cut from a lattice with
# infinite bonding, i.e. atomic solid, and (2) cut from a
# lattice with finite bonding, i.e. molecular or VDW solid, for
# (1) termination is relevant and is honored and so boundary
# information is necessary and molecules not completed, for
# (2) termination is not relevant and is not honored and so
# boundary information is not necessary and molecules are
# completed
if self.molecule_completion:
nanoparticle_indices = set(
range(1, self.nanoparticle.atom_total + 1))
while nanoparticle_indices:
aatom = self.nanoparticle.atom[nanoparticle_indices.pop()]
if not self.isInside(aatom):
continue
# make all atoms inside but label those that are outside
# so that the slab builder will treat them as if they
# are terminating atoms
for mol_atom in aatom.getMolecule().atom:
nanoparticle_indices.discard(mol_atom.index)
if not self.isInside(mol_atom):
mol_atom.property[ORIGINAL_ATOM_KEY] = False
self.outside_indices.remove(mol_atom.index)
self.inside_indices.add(mol_atom.index)
else:
for target_atom in self.nanoparticle.atom:
target_inside = self.isInside(target_atom)
boundary_indices = set()
for bonded_atom in target_atom.bonded_atoms:
bonded_inside = self.isInside(bonded_atom)
if bonded_inside != target_inside:
if self.term_frag != TERM_FRAGS[0]:
self.checkIntersectingBond(target_atom.index,
bonded_atom.index)
if bonded_inside:
boundary_index = (bonded_atom.index,
target_atom.index)
else:
boundary_index = (target_atom.index,
bonded_atom.index)
boundary_indices.add(boundary_index)
# update collections, if the target atom is inside and bonded exclusively
# to outside atoms then make the target atom an outside atom (this is
# a rare case but can happen for certain crystals when an inside atom
# is close to a polyhedron vertex)
if target_inside and target_atom.bond_total == len(
boundary_indices):
self.inside_indices.remove(target_atom.index)
self.outside_indices.add(target_atom.index)
self.boundary_indices.difference_update(boundary_indices)
else:
self.boundary_indices.update(boundary_indices)
# make unique lists
flattened = set()
for abond in self.boundary_indices:
flattened.update(set(abond))
self.inside_indices = sorted(self.inside_indices.difference(flattened))
self.outside_indices = sorted(
self.outside_indices.difference(flattened))
self.boundary_indices = list(self.boundary_indices)
[docs] def removeAndRenumber(self, to_remove=None, renumber_map=None):
"""
Remove the specified atoms from the nanoparticle and renumber all
collections of indices or just renumber the collections using the
specified renumber map.
:type to_remove: list or none
:param to_remove: indices of atoms to remove
:type renumber_map: dict or none
:param renumber_map: a renumber map to use in renumbering collections
"""
def process_list(alist):
return [renumber_map[idx] for idx in alist if \
renumber_map[idx] is not None]
# just return if there is nothing to do
if not to_remove and not renumber_map:
return
if to_remove:
renumber_map = self.nanoparticle.deleteAtoms(to_remove,
renumber_map=True)
boundary_indices_renumbered = []
for atuple in self.boundary_indices:
atuple_renumbered = tuple((renumber_map[x] for x in atuple))
if None not in atuple_renumbered:
boundary_indices_renumbered.append(atuple_renumbered)
self.boundary_indices = boundary_indices_renumbered
self.inside_indices = process_list(self.inside_indices)
self.outside_indices = process_list(self.outside_indices)
[docs] def cutAwayExcess(self):
"""
Remove the unneeded atoms from the crystalline supercell.
"""
# determine the atoms to remove, if there are boundary indices then
# handle splitting them up for two cases (1) no termination and (2)
# a specific type of termination in which certain atoms should not
# be terminated, if there are not boundary indices then just proceed
to_remove = list(self.outside_indices)
if self.boundary_indices:
if self.term_frag == TERM_FRAGS[0]:
keep, remove = list(zip(*(self.boundary_indices)))
self.inside_indices += list(set(keep))
to_remove += list(set(remove))
elif self.no_term_indices:
tmp_boundary_indices = list(self.boundary_indices)
self.boundary_indices = []
for atuple in tmp_boundary_indices:
inside, outside = atuple
if outside in self.no_term_indices:
if inside not in self.inside_indices:
self.inside_indices.append(inside)
if outside not in to_remove:
to_remove.append(outside)
else:
self.boundary_indices.append(atuple)
to_remove.sort()
# remove atoms and renumber collections
self.removeAndRenumber(to_remove=to_remove)
[docs] def handleOneToMany(self, inside_index, outside_indices):
"""
Handle the one to many situation, i.e. a single atom inside
the region of interest bonded to multiple atoms outside of
the region of interest.
:type inside_index: int
:param inside_index: the index of the inside atom
:type outside_indices: list
:param outside_indices: the indices of the outside atoms
"""
# remove bonding amongst outside atoms
to_remove = set()
for outside_index in outside_indices:
for outside_atom in self.nanoparticle.atom[
outside_index].bonded_atoms:
if outside_atom.index != inside_index:
to_remove.add(
tuple(sorted([outside_index, outside_atom.index])))
for atuple in to_remove:
self.nanoparticle.deleteBond(*atuple)
# make outside atoms hydrogen and set bond lengths
for index in outside_indices:
self.nanoparticle.atom[index].element = 'H'
self.nanoparticle.adjust(self.term_bond_length, inside_index, index)
self.nanoparticle.atom[index].property[ORIGINAL_ATOM_KEY] = False
[docs] def handleManyToOne(self, inside_indices, outside_index):
"""
Handle the many to one situation, i.e. multiple atoms
inside the region of interest are bonded to a single atom
that is outside of the region of interest.
:type inside_indices: list
:param inside_indices: the indices of the inside atoms
:type outside_index: int
:param outside_index: the index of the outside atom
:rtype: list
:return: contains (inside, outside) tuples of newly created
bonds
"""
outside_vec = numpy.array(self.nanoparticle.atom[outside_index].xyz)
# terminate each inside index by bonding to a newly added hydrogen
# along the pre-existing bond vectors, clashes will be resolved with
# subsequent optimization of the structure with a real physical model
new_bonds = []
for inside_index in inside_indices:
inside_vec = numpy.array(self.nanoparticle.atom[inside_index].xyz)
bond_unit_vec = transform.get_normalized_vector(outside_vec -
inside_vec)
x_term, y_term, z_term = inside_vec + self.term_bond_length * bond_unit_vec
h_term = self.nanoparticle.addAtom('H', x_term, y_term, z_term)
self.nanoparticle.deleteBond(inside_index, outside_index)
self.nanoparticle.addBond(inside_index, h_term.index, 1)
new_bonds.append((inside_index, h_term.index))
return new_bonds
[docs] def configureBonds(self):
"""
Configure the nanoparticle bonds that intersect the shape.
"""
# update data structure to ease collecting of intersecting bonds
# that share an inside or outside atom
boundary_indices_grouped = \
group_subtuples_by_key_and_value(tuple(self.boundary_indices))
# handle three situations for nanoparticles: (1) one-to-many, a single
# atom in the region of interest is bonded to multiple atoms outside
# the region of interest, here the outside atoms themselves are used
# to terminate, (2) many-to-one, multiple atoms in the region of
# interest are bonded to a single atom outside the region of
# interest, here the outside atom can't be used and so we create
# several new terminating atoms, (3) one-to-one, simple case
to_remove = []
for atuple in boundary_indices_grouped:
inside, outside = atuple
if len(inside) > 1:
new_bonds = self.handleManyToOne(list(inside), outside[0])
obsolete_bonds = set(zip(inside, len(inside) * outside))
old_bonds = list(
set(self.boundary_indices).difference(obsolete_bonds))
self.boundary_indices = old_bonds + new_bonds
to_remove.append(outside[0])
else:
self.handleOneToMany(inside[0], list(outside))
# remove atoms and renumber collections
self.removeAndRenumber(to_remove=to_remove)
[docs] def doTermination(self):
"""
Terminate the nanoparticle.
"""
def get_bonded_indices(aatom):
return [bonded_atom.index for bonded_atom in \
self.nanoparticle.atom[aatom].bonded_atoms]
frag_group, frag_name = TERM_DICT[self.term_frag]
term_bonds = list(self.boundary_indices)
while term_bonds:
# get from and to indices
from_index, to_index = term_bonds.pop()
# get the other bonds for the from atom (this is so that we can later
# figure out the atom index of the fragment atom that has replaced
# the terminating hydrogen)
bonded_atoms = get_bonded_indices(from_index)
bonded_atoms.remove(to_index)
# attach the fragment
renumber_map = build.attach_fragment(self.nanoparticle, from_index,
to_index, frag_group,
frag_name, TERM_DIRECTION)
# determine the atom number of the fragment atom that has replaced
# the terminating hydrogen
bonded_atoms = set(
[renumber_map[bonded_atom] for bonded_atom in bonded_atoms])
updated_to_index = \
set(get_bonded_indices(renumber_map[from_index])).difference(bonded_atoms).pop()
# adjust the bond length of the new nanoparticle-fragment bond
self.nanoparticle.adjust(self.term_bond_length, renumber_map[from_index], \
updated_to_index)
# update term_bonds
term_bonds = [
(renumber_map[x[0]], renumber_map[x[1]]) for x in term_bonds
]
# update collections but first add in the new terminating bond in terms
# of the original indexing, i.e. the keys of the renumber map
self.boundary_indices.append(
(from_index, self.nanoparticle.atom_total + 1))
self.removeAndRenumber(renumber_map=renumber_map)
[docs] def setNanoparticleProperties(self):
"""
Set nanoparticle properties.
"""
# set title
self.nanoparticle.title = self.supercell.title + '-' + self.shape
# set atom and bond visualization representation
set_representation(self.nanoparticle)
# set number of nanoparticle atoms
if self.boundary_indices:
inside_indices = list(set(list(zip(*(self.boundary_indices)))[0]))
else:
inside_indices = []
inside_indices = self.inside_indices + inside_indices
self.nanoparticle_natoms = len(inside_indices)
self.nanoparticle.property[
NANOPARTICLE_NATOMS_KEY] = self.nanoparticle_natoms
# set nanoparticle volume
self.nanoparticle_volume = self.shape_object.volume
self.nanoparticle.property[
NANOPARTICLE_VOLUME_KEY] = self.nanoparticle_volume
# set nanoparticle surface area
self.nanoparticle_surface_area = self.shape_object.surface_area
self.nanoparticle.property[NANOPARTICLE_SURFACE_AREA_KEY] = \
self.nanoparticle_surface_area
# set nanoparticle circumradius
if self.polyhedron and self.shape_object.TYPE != shapes.PRISM:
self.nanoparticle_circumradius = self.shape_object.circumradius
self.nanoparticle.property[NANOPARTICLE_CIRCUMRADIUS_KEY] = \
self.nanoparticle_circumradius
# set nanoparticle density
volume_cm3 = self.nanoparticle_volume * math.pow(
old_div(1.0, CM_TO_ANGSTROM), 3)
total_mass = 0.0
for index in inside_indices:
total_mass += self.nanoparticle.atom[index].atomic_weight
total_mass = old_div(total_mass, scipy_constants.N_A)
self.nanoparticle_density = old_div(total_mass, volume_cm3)
self.nanoparticle.property[NANOPARTICLE_DENSITY_KEY] = \
self.nanoparticle_density
[docs] def printParams(self):
"""
Log the parameters.
"""
HEADER = 'Nanoparticle Parameters'
# shape parameters
shape_name = textlogger.get_param_string('Shape name', self.shape,
self.MSGWIDTH)
shape_type = textlogger.get_param_string('Shape type',
self.shape_object.TYPE,
self.MSGWIDTH)
num_unique_faces = textlogger.get_param_string(
'Num. primary axes', self.shape_object.NUM_UNIQUE_FACES,
self.MSGWIDTH)
num_unique_edges = textlogger.get_param_string(
'Num. secondary axes', self.shape_object.NUM_UNIQUE_EDGES,
self.MSGWIDTH)
shape_origin = textlogger.get_param_string('Shape origin', self.origin,
self.MSGWIDTH)
primary_alignment = textlogger.get_param_string(
'Primary alignment', dict_to_str(self.primary_alignment_dict),
self.MSGWIDTH)
secondary_alignment = textlogger.get_param_string(
'Secondary alignment', dict_to_str(self.secondary_alignment_dict),
self.MSGWIDTH)
overwrite_vertices = textlogger.get_param_string(
'Overwrite vertices given', bool(len(self.overwrite_vertices)),
self.MSGWIDTH)
shape_params = [
'%s (%s)' % (str(x[0]), x[1])
for x in zip(self.params, self.shape_object.DESCRIPTION)
]
shape_params = ', '.join(shape_params)
shape_params = textlogger.get_param_string(
'Shape parameters / Ang. & degree', shape_params, self.MSGWIDTH)
shape_volume = '%.4f' % round(self.nanoparticle_volume, NUM_DECIMAL)
shape_volume = textlogger.get_param_string('Shape volume / Ang.^3',
shape_volume, self.MSGWIDTH)
shape_surface_area = '%.4f' % round(self.nanoparticle_surface_area,
NUM_DECIMAL)
shape_surface_area = textlogger.get_param_string(
'Shape surface area / Ang.^2', shape_surface_area, self.MSGWIDTH)
if self.polyhedron:
shape_num_verts = textlogger.get_param_string(
'Shape number of vertices', len(self.shape_object.vertices),
self.MSGWIDTH)
shape_num_faces = textlogger.get_param_string(
'Shape number of faces', len(self.shape_object.faces),
self.MSGWIDTH)
if self.shape_object.TYPE != shapes.PRISM:
shape_circumradius = '%.4f' % round(
self.nanoparticle_circumradius, NUM_DECIMAL)
shape_circumradius = textlogger.get_param_string(
'Shape circumradius / Ang.', shape_circumradius,
self.MSGWIDTH)
embedded = textlogger.get_param_string('Embedded in supercell', \
self.embedded, self.MSGWIDTH)
# termination parameters
term_frag = textlogger.get_param_string('Terminating fragment',
self.term_frag, self.MSGWIDTH)
term_bond_length = textlogger.get_param_string(
'Terminating bond length / Ang.', self.term_bond_length,
self.MSGWIDTH)
allow_fragments = textlogger.get_param_string('Allow fragments',
self.allow_fragments,
self.MSGWIDTH)
# nanoparticle parameters
nanoparticle_natoms = textlogger.get_param_string(
'Nanoparticle number of atoms', self.nanoparticle_natoms,
self.MSGWIDTH)
nanoparticle_density = '%.4f' % round(self.nanoparticle_density,
NUM_DECIMAL)
nanoparticle_density = textlogger.get_param_string(
'Nanoparticle density / g/cm^3', nanoparticle_density,
self.MSGWIDTH)
# lattice parameters
if self.lattice_parameters:
lattice_parameters = \
textlogger.get_param_string('Lattice (a,b,c,alpha,beta,gamma) / Ang. & degree',
', '.join([str(round(param, 3)) for param in self.lattice_parameters]),
self.MSGWIDTH)
self.logger.info(HEADER)
self.logger.info('-' * len(HEADER))
self.logger.info(shape_name)
self.logger.info(shape_type)
self.logger.info(num_unique_faces)
self.logger.info(num_unique_edges)
self.logger.info(shape_origin)
self.logger.info(primary_alignment)
self.logger.info(secondary_alignment)
self.logger.info(overwrite_vertices)
self.logger.info(shape_params)
self.logger.info(shape_volume)
self.logger.info(shape_surface_area)
if self.polyhedron:
self.logger.info(shape_num_verts)
self.logger.info(shape_num_faces)
if self.shape_object.TYPE != shapes.PRISM:
self.logger.info(shape_circumradius)
self.logger.info(embedded)
self.logger.info(term_frag)
self.logger.info(term_bond_length)
self.logger.info(allow_fragments)
self.logger.info(nanoparticle_natoms)
self.logger.info(nanoparticle_density)
if self.lattice_parameters:
self.logger.info(lattice_parameters)
self.logger.info('')
[docs] def getAlignmentVector(self, alignment_dict):
"""
Return the alignment vector.
:type alignment_dict: dict
:param alignment_dict: defines the alignment
:rtype: numpy.array
:return: the vector on which to align
"""
NO_KEYS_MSG = (
'You have specified an alignment using an %s or %s basis. '
'Such alignments require the following structure properties to exist '
'in your input Maestro file of the supercell, i.e. %s.')
if alignment_dict[BASIS_KEY] != XYZ_VALUE:
lattice_parameter_keys = get_lattice_properties()
if not check_if_lattice_properties_exist(
self.supercell, keys=lattice_parameter_keys):
msg = NO_KEYS_MSG % (HKL_VALUE, ABC_VALUE,
', '.join(lattice_parameter_keys))
if self.logger:
self.logger.error(msg)
raise ValueError(msg)
else:
self.lattice_parameters = [
self.supercell.property[x] for x in lattice_parameter_keys
]
a_vec, b_vec, c_vec = xtal.get_lattice_vectors(
*(self.lattice_parameters))
if alignment_dict[BASIS_KEY] == ABC_VALUE:
terms = list(
zip(alignment_dict[VECTOR_KEY], [a_vec, b_vec, c_vec]))
alignment_vector = sum(
[coef * vector for coef, vector in terms])
elif alignment_dict[BASIS_KEY] == HKL_VALUE:
h_index, k_index, l_index = alignment_dict[VECTOR_KEY]
self.plane = plane.CrystalPlane(h_index, k_index, l_index,
a_vec, b_vec, c_vec)
alignment_vector = self.plane.normal_vec
else:
alignment_vector = alignment_dict[VECTOR_KEY]
return alignment_vector
[docs] def translateToPlane(self):
"""
Translate the shape object so that the reference face is
coplanar with the requested crystal plane.
"""
face_idx = self.primary_alignment_dict[AXIS_KEY]
face = self.shape_object.reference_faces[face_idx - 1]
unit_normal = transform.get_normalized_vector(self.plane.normal_vec)
to_plane = self.plane.normal_vec - numpy.dot(face.center,
unit_normal) * unit_normal
self.shape_object.translateVertices(self.shape_object.vertices,
to_plane)
[docs] def overwriteVertices(self):
"""
Overwrite the vertices of the shape object.
"""
BAD_VERTS_MSG = (
'You must specify at least 4 overwrite vertices and each '
'vertex must be defined using 3 floats.')
if len(self.overwrite_vertices) < 12 or len(
self.overwrite_vertices) % 3:
if self.logger:
self.logger.error(BAD_VERTS_MSG)
raise ValueError(BAD_VERTS_MSG)
vertices = []
for index, first_index in enumerate(
list(range(0, len(self.overwrite_vertices), 3)), 1):
coord = numpy.array(
self.overwrite_vertices[first_index:first_index + 3])
vertices.append(shapes.Vertex(index, coord))
self.shape_object.updateShape(vertices)
[docs] def runIt(self, template_only=False):
"""
Create the nanoparticle.
:type template_only: bool
:param template_only: If True, the method will return immediately after
the template is created and before the nanoparticle has been fully
synthesized. Using this parameter means that the .template structure
is valid but no .nanoparticle structure is created - the .nanoparticle
property will be set to None.
"""
# check input
self.checkInput()
# get reference face and edge axis IDs
ref_face_idx = self.primary_alignment_dict[AXIS_KEY]
ref_edge_idx = self.secondary_alignment_dict[AXIS_KEY]
# get along vectors
ref_face_normal_along = self.getAlignmentVector(
self.primary_alignment_dict)
ref_edge_along = self.getAlignmentVector(self.secondary_alignment_dict)
# get shape class
shape_object = shapes.get_shape_object_by_name(self.shape)
# set whether or not we have a polyhedral shape
self.polyhedron = shape_object.TYPE in shapes.POLYHEDRON_TYPES
# get shape object
try:
if self.polyhedron:
self.shape_object = shape_object(*(self.params), center=self.origin, \
ref_face_idx=ref_face_idx, ref_face_normal_along=ref_face_normal_along, \
ref_edge_idx=ref_edge_idx, ref_edge_along=ref_edge_along)
else:
self.shape_object = shape_object(*(self.params),
center=self.origin)
except ValueError as msg:
if self.logger:
self.logger.error(msg)
raise msg
# if we have a polyhedron then do some things
if self.polyhedron:
# translate reference face to crystal plane along the normal and allow the
# user specified origin to change in this case
if self.primary_alignment_dict[BASIS_KEY] == HKL_VALUE:
self.translateToPlane()
self.origin = self.shape_object.center
# overwrite vertices
if self.overwrite_vertices:
self.overwriteVertices()
# set template
if self.include_axes:
self.shape_object.addAlignmentAxesToTemplate()
self.template = self.shape_object.template.copy()
# Mark the template
if self.template:
self.template.property[NANOPARTICLE_TEMPLATE] = True
if template_only:
# The nanoparticle isn't valid at this point, so blank it out
self.nanoparticle = None
return
# determine how to handle molecules at the nanoparticle boundary
self.setMoleculeCompletion()
xtal.delete_all_pbc_bonds(self.nanoparticle)
# the following issue pertains to convex polyhedron shapes only,
# in shapes.Face.intersectSegmentAndPlane, from the following execution
# path partitionAtoms -> isInside -> shapes.<Shape>.pointInside, a
# divide-by-zero RuntimeWarning can be shown, we want to silent that
# because it is being properly processed in that intersectSegmentAndPlane
# method
warnings.simplefilter('ignore', RuntimeWarning)
# partition atoms
self.partitionAtoms()
# check if we are completely outside the region of interest
if len(self.outside_indices) == self.nanoparticle.atom_total:
NO_ATOMS_MSG = (
'Given your provided crystalline supercell '
'and shape parameters this nanoparticle contains zero '
'atoms.')
if self.logger:
self.logger.error(NO_ATOMS_MSG)
raise ValueError(NO_ATOMS_MSG)
# if we have a polyhedron then check if it was embedded in the
# supercell such that all faces, not planes containing faces, have been
# intersected by at least one pointInside query, this information
# tells us if part of the requested polyhedron lies outside the
# region of interest, this is currently allowed as it allows
# users to be more flexible with the way in which they use this
# script, for example they can pit surfaces, etc.
if self.polyhedron:
self.embedded = \
self.shape_object.allFacesIntersected()
if not self.embedded:
if self.logger:
NOT_EMBEDDED_MSG = (
'The polyhedron that has been '
'specified is not entirely embedded in the '
'given crystalline supercell, meaning that '
'there is at least one face hanging outside '
'of the supercell boundaries. The returned '
'nanoparticle may thus be incomplete, i.e. '
'will not entirely fill the shape you requested.')
self.logger.warning(NOT_EMBEDDED_MSG)
# cut away unneeded atoms and if terminating then configure
# the bonds and terminate the nanoparticle
self.cutAwayExcess()
if self.term_frag != TERM_FRAGS[0]:
self.configureBonds()
if self.term_frag != TERM_FRAG:
self.doTermination()
color.apply_color_scheme(self.nanoparticle, 'element')
self.nanoparticle.retype()
# set nanoparticle properties
self.setNanoparticleProperties()
# log nanoparticle properties
if self.logger:
self.printParams()
[docs]def build_kwargs(kwargs_list,
defaults,
acceptable_keys=None,
token_separator=TOKEN_SEPARATOR,
key_value_separator=KEY_VALUE_SEPARATOR,
all_choices=None):
"""
Build the kwargs dictionary from an unformatted and unvalidated
list of kwargs.
:type kwargs_list: list
:param kwargs_list: list of strings that contain kwargs in
one form or another
:type defaults: dict
:param defaults: default kwargs to be updated
:type acceptable_keys: list
:param acceptable_keys: list of acceptable keys
:type token_separator: str
:param token_separator: the separator used for tokens
:type key_value_separator: str
:param key_value_separator: the separator used in key-value pairs
:type all_choices: dict
:param all_choices: if the values of certain keys are limited then
collect those options here in lists for validation purposes
:rtype: dict
:return: dictionary of kwargs, formatted, validated, and contains
unspecified defaults
"""
if acceptable_keys is None:
acceptable_keys = ALIGNMENT_KEYS
if all_choices is None:
all_choices = {BASIS_KEY: BASIS_VALUES}
general_msg = ('Alignment options must be specified using tokens of the form '
'<key>%s<value>, where key is a valid key from %s, and where multiple '
'tokens are separated with a \'%s\'.') % \
(key_value_separator, acceptable_keys, token_separator)
choices_msg = ('The %s sub-option must take on a value in %s.')
pairs = []
for token in ' '.join(kwargs_list).split(token_separator):
if token.count(key_value_separator):
key, value = [
x.strip() for x in token.strip().split(key_value_separator, 1)
]
if key not in acceptable_keys or value == '':
raise ValueError(general_msg)
else:
choices = all_choices.get(key)
if choices is None or value in choices:
pairs.append((key, value))
else:
raise ValueError(choices_msg % (key, choices))
else:
raise ValueError(general_msg)
kwargs = defaults.copy()
kwargs.update(OrderedDict(pairs))
return kwargs
[docs]def set_representation(astructure):
"""
Set the representation of the given structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure to set the representation for
"""
astructure.applyStyle()