"""
Classes and functions for making honeycomb nanosheets.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
import math
from collections import OrderedDict
from past.utils import old_div
import numpy
import schrodinger.structure as structure
from schrodinger.application.matsci import textlogger as tlog
from schrodinger.application.matsci.coarsegrain import set_atom_coarse_grain_properties
from schrodinger.application.matsci.coarsegrain import set_as_coarse_grain
from schrodinger.application.matsci.nano import check
from schrodinger.application.matsci.nano import constants
from schrodinger.application.matsci.nano import slab
from schrodinger.application.matsci.nano import util
from schrodinger.application.matsci.nano import xtal
from schrodinger.structutils import build
from schrodinger.structutils import color
from schrodinger.structutils import minimize
from schrodinger.structutils import transform
from schrodinger.test import ioredirect
_version = '$Revision 0.0 $'
[docs]class SymmetryEquiv(object):
"""
Manage symmetry equivalent positions in the unit cell.
"""
[docs] def __init__(self, element, genvec, indicies, is_coarse_grain=False):
"""
Create an instance.
:type element: str
:param element: atom for this symmetry type
:type genvec: numpy.array
:param genvec: vector to generate symmetry equivalent positions
:type indicies: list of ints
:param indicies: atom indicies to use for the symmetry equivalent
atoms
:type is_coarse_grain: bool
:param is_coarse_grain: Whether a coarse grain structure is being
created
"""
self.is_coarse_grain = is_coarse_grain
self.element = (constants.Constants.CG_AA_ELEMENT_SYMBOL
if self.is_coarse_grain else element)
self.name = element if self.is_coarse_grain else None
self.genvec = genvec
self.indicies = list(indicies)
self.positions = {self.indicies[0]: self.genvec}
self.generateSymmetryEquiv()
[docs] def generateSymmetryEquiv(self):
"""
Generate symmetry equivalent positions in the unit cell.
"""
genvec = self.genvec
for index in self.indicies[1:]:
genvec = util.get_rotated_vector(genvec, \
HoneycombUnitCell.ANGLELARGE)
self.positions[index] = genvec
[docs]class HoneycombUnitCell(object):
r"""
Create a unit cell of a honeycomb lattice. The unit cell will be a
regular hexagon according to the following coordinate system::
-----------------------------------------
| Y |
| -------------------- |
| / 6 5 \ |
| / \ |
| / \ |
| / \ |
| / \ |
| / \ |
| | 1 O 4 | X |
| \ / |
| \ / |
| \ / |
| \ / |
| \ / |
| \ 2 3 / |
| -------------------- |
-----------------------------------------
The hexagon will be centered on the origin O where vec(X, O) and vec(Y, O)
are the X- and Y-axes and the Z-coordinate is zero. All edge lengths are
equivalent and internal small triangles, i.e. triangle(1, O, 2) are
equilateral. The following angles will be important: angle(1, 2, 3) =
120 degrees, angle(1, 2, O) = 60 degrees, and angle(1, 2, 6) = 30 degrees.
In general the asymmetric unit will contain two atoms. The symmetry
operation is C3 and thus symmetry equivalent positions fall into two sets,
i.e. [1, 3, 5] and [2, 4, 6]. This unit cell will have single bonds assigned.
"""
ANGLELARGE = math.radians(120)
ANGLEMEDIUM = math.radians(60)
ANGLESMALL = math.radians(30)
NUMATOMS = 6
NUM_EQUIV_SETS = 2
NUM_SYM_ATOMS = old_div(NUMATOMS, NUM_EQUIV_SETS)
SYMATOMS = [[1, 3, 5], [2, 4, 6]]
UNITCELLINDEX = 'i_matsci_Unit_Cell_Index'
[docs] def __init__(self, element1, element2, bondlength, is_coarse_grain=False):
"""
Create an instance.
:type element1: str
:param element1: elemental symbol of the first atom
:type element2: str
:param element2: elemental symbol of the second atom
:type bondlength: float
:param bondlength: bond length between the first and second
atoms in Angstrom
:type is_coarse_grain: bool
:param is_coarse_grain: Whether a coarse grain structure is being
created
"""
self.element1 = element1
self.element2 = element2
self.bondlength = bondlength
self.is_coarse_grain = is_coarse_grain
self.structure = None
self.atomic_number1 = None
self.atomic_number2 = None
[docs] def getAsymmetricUnit(self):
"""
Return two vectors to the atoms in the asymmetric unit.
:rtype: numpy.array, numpy.array
:return: asupos1, asupos2, positions of the two atoms in the
asymmetric unit
"""
asupos1 = numpy.array([-1 * self.bondlength, 0, 0])
asupos2 = util.get_rotated_vector(asupos1, \
self.ANGLEMEDIUM)
return asupos1, asupos2
[docs] def getStructure(self, symgroup1, symgroup2):
"""
Create Structure for the unit cell.
:type symgroup1: SymmetryEquiv
:param symgroup1: first equivalent group
:type symgroup2: SymmetryEquiv
:param symgroup2: second equivalent group
:rtype: schrodinger.structure.Structure
:return: unitcell, structure for the unit cell
"""
# create structure by adding symmetry equivalent sets of atoms,
# create element and xyz properties, setting the element here
# can throw warnings (see MATSCI-2880) so silence it
unit_cell_structure = structure.create_new_structure(self.NUMATOMS)
for symset in [symgroup1, symgroup2]:
for index, position in symset.positions.items():
current_atom = unit_cell_structure.atom[index]
with ioredirect.IOSilence():
current_atom.element = symset.element
current_atom.xyz = position
if symset.is_coarse_grain:
# Assume covalent radius is same as the radius of the
# coarse grain particle.
set_atom_coarse_grain_properties(unit_cell_structure,
current_atom,
symset.name,
radius=self.bondlength / 2)
# create bonds, make sure to close the ring
for index in range(2, self.NUMATOMS + 1):
atoma = unit_cell_structure.atom[index - 1]
atomb = unit_cell_structure.atom[index]
unit_cell_structure.addBond(atoma, atomb, 1)
unit_cell_structure.addBond(unit_cell_structure.atom[1], atomb, 1)
# do atom types and color by element
unit_cell_structure.retype()
color.apply_color_scheme(unit_cell_structure, 'element')
return unit_cell_structure
[docs] def setAtomProps(self):
"""
Create structure.atom properties for later use.
"""
# this properties is used in connecting new unit cells
# to the lattice while growing it. This property will
# be removed from the structure.atoms once the lattice
# has been formed. This property gives the unit cell
# numbering.
for atom in self.structure.atom:
atom.property[self.UNITCELLINDEX] = atom.index
[docs] def buildUnitCell(self):
"""
Build the unit cell.
"""
asuvec1, asuvec2 = self.getAsymmetricUnit()
symgroup1 = SymmetryEquiv(self.element1,
asuvec1,
self.SYMATOMS[0],
is_coarse_grain=self.is_coarse_grain)
symgroup2 = SymmetryEquiv(self.element2,
asuvec2,
self.SYMATOMS[1],
is_coarse_grain=self.is_coarse_grain)
self.structure = self.getStructure(symgroup1, symgroup2)
self.setAtomProps()
self.atomic_number1, self.atomic_number2 = \
[self.structure.atom[x].atomic_number for x in [1, 2]]
[docs]class NeighborType(object):
"""
Manage the properties of a neighbor type.
"""
[docs] def __init__(self, edgetype, bondingatoms):
"""
Create an instance.
:type edgetype: str
:param edgetype: neighbor connection type
:type bondingatoms: list of ints
:param bondingatoms: bonding atom indices, in unit cell numbering,
of the neighbor
"""
self.edgetype = edgetype
self.bondingatoms = list(bondingatoms)
[docs]class Grow(object):
"""
Manage the details of how to attach cells to the lattice of the desired
shape. Most of this class is really just to manage grow parameters
which are derived by applying the symmetry operations on a set of basic
grow parameters given below.
"""
# here are three sets in unit cell numbering (1) terminating indicies,
# (2) growing indicies, and (3) all indicies. See from the HoneycombUnitCell
# class that terminating indicies are terminating (outwards facing, i.e.
# marking the end of the lattice) for the first cell placed in the lattice,
# i.e. the left corner. And the inverse of which is terminating for the
# last cell placed in the lattice, i.e. the right corner. The opposite
# is true of the grow indicies.
TERMINDICIES = [1, 2, 6]
GROWINDICIES = [util.get_inversion_index(index) for index in \
TERMINDICIES]
GROWINDICIES.sort()
ALLINDICIES = TERMINDICIES + GROWINDICIES
ALLINDICIES.sort()
# along the lattice edges we have the following base sets of terminating
# indicies
ZIGZAG_BY_ONE_BASE = [2]
ARMCHAIR_BY_ONE_BASE = [ZIGZAG_BY_ONE_BASE[0] - 1, ZIGZAG_BY_ONE_BASE[0]]
ONE_BY_ZIGZAG_BASE = [6]
ONE_BY_ARMCHAIR_BASE = [1, ONE_BY_ZIGZAG_BASE[0]]
# the lattice is grown as a matrix row-wise and so a given cell can
# be indexed by two integers as multiples of the two grow vectors,
# a given cell can then potentially has neighbors at the following
# positions relative to the target cell. Neighbors share either
# an edge or an atom. Along a zigzag direction each cell placed has
# been explicity created (edge-to-edge), while along an armchair direction
# cells are only placed (atom-to-atom), and so many of the cells in the
# lattice are actually formed simply out of the empty space left behind.
TYPE1NEIGHBOR = (-1, 0)
TYPE2NEIGHBOR = (0, -1)
TYPEZZNEIGHBOR = (-1, 1)
TYPEAANEIGHBOR = (-1, -1)
[docs] def __init__(self, edgetype1, edgetype2, ncell1, ncell2):
"""
Create an instance.
:type edgetype1: str
:param edgetype1: type of edge for lattice side 1
:type edgetype2: str
:param edgetype2: type of edge for lattice side 2
:type ncell1: int
:param ncell1: number of cells along lattice side 1
:type ncell2: int
:param ncell2: number of cells along lattice side 2
"""
self.edgetype1 = edgetype1
self.edgetype2 = edgetype2
self.ncell1 = ncell1
self.ncell2 = ncell2
[docs] def getGrowVectors(self, lattvec1, lattvec2):
"""
Return the two vectors needed to grow the honeycomb lattice with
the specified types of edges.
:type lattvec1: numpy.array
:param lattvec1: lattice vector 1
:type lattvec2: numpy.array
:param lattvec2: lattice vector 2
:rtype: numpy.array, numpy.array
:return: growvec1, growvec2, the two lattice grow vectors
"""
# the lattice vectors form a basis for the grow vectors
# which are the vectors which are actually used to grow the
# lattice of the desired shape
if self.edgetype1 == constants.Constants.ARMCHAIR:
growvec1 = 2 * lattvec1 - lattvec2
else:
growvec1 = lattvec1
if self.edgetype2 == constants.Constants.ARMCHAIR:
growvec2 = 2 * lattvec2 - lattvec1
else:
growvec2 = lattvec2
return growvec1, growvec2
[docs] def getNeighborTypes(self):
"""
Return a dictionary of neighbor types.
:rtype: dict
:return: neighbortypes, as keys has tuples of deltas to
neighboring positions and as values has NeighborTypes
"""
# create the dictionary which maps potential neighboring
# cells to the type of edge connecting them to the target
# cell as well as the the list of atoms which would be
# bonded to the target cell
neighbortypes = {}
delta = self.TYPE1NEIGHBOR
bondingatoms = self.GROWINDICIES[:-1]
if self.edgetype1 == constants.Constants.ARMCHAIR:
bondingatoms = bondingatoms[:-1]
neighbortype = NeighborType(self.edgetype1, bondingatoms)
neighbortypes[delta] = neighbortype
delta = self.TYPE2NEIGHBOR
bondingatoms = self.GROWINDICIES[1:]
if self.edgetype2 == constants.Constants.ARMCHAIR:
bondingatoms = bondingatoms[1:]
neighbortype = NeighborType(self.edgetype2, bondingatoms)
neighbortypes[delta] = neighbortype
if self.edgetype1 == self.edgetype2:
if self.edgetype1 == constants.Constants.ZIGZAG:
delta = self.TYPEZZNEIGHBOR
bondingatoms = [2, self.GROWINDICIES[0]]
else:
delta = self.TYPEAANEIGHBOR
bondingatoms = [self.GROWINDICIES[1]]
neighbortype = NeighborType(self.edgetype1, bondingatoms)
neighbortypes[delta] = neighbortype
return neighbortypes
[docs] def getRedundantAtoms(self, deltas, neighbortypes):
"""
Return a list of atoms in the current cell that will be redundant
with its neighboring cells when placed in the lattice. Unit cell
numbering, rather than lattice numbering, is used.
:type deltas: list of tuples
:param deltas: deltas to neighboring cells
:type neighbortypes: dict
:param neighbortypes: contains as keys the deltas for potential
neighbors and as values NeighborTypes
:rtype: list of ints
:return: redundantatoms, list of redundant atoms
"""
redundantatoms = set()
for delta in deltas:
neighbortype = neighbortypes[delta]
if neighbortype.edgetype != constants.Constants.ARMCHAIR:
for bondingatom in neighbortype.bondingatoms:
redundantatom = \
util.get_inversion_index(bondingatom)
redundantatoms.add(redundantatom)
redundantatoms = list(redundantatoms)
redundantatoms.sort()
return redundantatoms
[docs] def getBondingInfo(self, deltas, neighbortypes):
"""
Return the details for attaching cells to the lattice.
:type deltas: list of tuples
:param deltas: deltas to neighboring cells
:type neighbortypes: dict
:param neighbortypes: contains as keys the deltas for potential
neighbors and as values NeighborTypes
:rtype: dict
:return: bonds_to_neighbors, dictionary containing as keys
the deltas to neighboring cells from cell (growindex1, growindex2)
and as values a list of neighbor-cell atom pairs to bond. Pairs
are given using the unit cell numbering rather than the lattice
numbering.
"""
bonds_to_neighbors = {}
# loop over neighbors, find neighbor-cell bonding atom pairs, make
# bonds handling special cases
for delta in deltas:
neighbortype = neighbortypes[delta]
cell_bonding_atoms = []
for bondingatom in neighbortype.bondingatoms:
invatom = util.get_inversion_index(bondingatom)
cell_bonding_atoms.append(invatom)
if neighbortype.edgetype == constants.Constants.ZIGZAG:
potential_bonding_atoms = \
set(self.ALLINDICIES).difference(set(cell_bonding_atoms))
potential_bonding_atoms = list(potential_bonding_atoms)
cell_bonding_atoms = \
[min(potential_bonding_atoms), max(potential_bonding_atoms)]
bonds = []
for natom, catom in zip(neighbortype.bondingatoms,
cell_bonding_atoms):
bond = [natom, catom]
bonds.append(bond)
bonds_to_neighbors[delta] = bonds
return bonds_to_neighbors
[docs] def getAtomsToTerminate(self, growindex1, growindex2, logger=None):
"""
Return a list of atoms in the current cell that are not shared with
any other cell. Unit cell numbering, rather than lattice numbering,
is used.
:type growindex1: int
:param growindex1: grow index 1
:type growindex2: int
:param growindex2: grow index 2
:type logger: logging.getLogger
:param logger: output logger
:rtype: list of ints
:return: atoms_to_terminate, list of unshared atoms
"""
def get_linear_termination(basetype1, basetype2, edgetype, ncells,
growindex):
if edgetype == constants.Constants.ARMCHAIR:
indicies = set(basetype2)
else:
indicies = set(basetype1)
invindicies = set([util.get_inversion_index(index) \
for index in indicies])
indicies.update(invindicies)
if not growindex:
indicies = indicies.union(set(self.TERMINDICIES))
elif growindex == ncells - 1:
indicies = indicies.union(set(self.GROWINDICIES))
indicies = list(indicies)
indicies.sort()
return indicies
def get_corner_termination(growindex):
if self.edgetype1 == constants.Constants.ARMCHAIR:
indicies1 = list(self.ARMCHAIR_BY_ONE_BASE)
else:
indicies1 = list(self.ZIGZAG_BY_ONE_BASE)
if self.edgetype2 == constants.Constants.ARMCHAIR:
indicies2 = list(self.ONE_BY_ARMCHAIR_BASE)
else:
indicies2 = list(self.ONE_BY_ZIGZAG_BASE)
if growindex:
indicies2 = [util.get_inversion_index(index) \
for index in indicies2]
else:
indicies1 = [util.get_inversion_index(index) \
for index in indicies1]
indicies1.extend(indicies2)
indicies1.sort()
return indicies1
def get_edge_termination(basetype1, basetype2, edgetype, ncells,
growindex):
if edgetype == constants.Constants.ARMCHAIR:
indicies = list(basetype2)
else:
indicies = list(basetype1)
if growindex == ncells - 1:
indicies = [util.get_inversion_index(index) \
for index in indicies]
indicies.sort()
return indicies
atoms_to_terminate = []
# the 1-by-1 case is supported mainly in case a user wants to
# make bilayers out of it, 1D lattices 1-by-N and N-by-1 are also
# supported, as well as obviously the 2D lattice N-by-N, in which
# case we have corner- and edge-type terminations
# single unit cell
if self.ncell1 == 1 and self.ncell2 == 1:
atoms_to_terminate = list(self.ALLINDICIES)
# 1D lattice along lattice dimension 1
elif self.ncell1 != 1 and self.ncell2 == 1:
atoms_to_terminate = get_linear_termination(
self.ZIGZAG_BY_ONE_BASE, self.ARMCHAIR_BY_ONE_BASE,
self.edgetype1, self.ncell1, growindex1)
# 1D lattice along lattice dimension 2
elif self.ncell1 == 1 and self.ncell2 != 1:
atoms_to_terminate = get_linear_termination(
self.ONE_BY_ZIGZAG_BASE, self.ONE_BY_ARMCHAIR_BASE,
self.edgetype2, self.ncell2, growindex2)
# 2D lattice
elif self.ncell1 != 1 and self.ncell2 != 1:
# do left, right, top, and bottom corners
# left corner, corner closest to the origin, beginning growth
if growindex1 == 0 and growindex2 == 0:
atoms_to_terminate = list(self.TERMINDICIES)
# right corner, corner opposite the origin, ending growth
elif growindex1 == self.ncell1 - 1 and growindex2 == self.ncell2 - 1:
atoms_to_terminate = list(self.GROWINDICIES)
# top and bottom corners respectively
elif growindex1 == 0 and growindex2 == self.ncell2 - 1 or \
growindex1 == self.ncell1 - 1 and growindex2 == 0:
atoms_to_terminate = get_corner_termination(growindex1)
# do left, right, top, and bottom lattice sides
# top and bottom lattice sides
elif growindex1 == 0 or growindex1 == self.ncell1 - 1 and \
growindex2 != self.ncell2 - 1:
atoms_to_terminate = get_edge_termination(
self.ONE_BY_ZIGZAG_BASE, self.ONE_BY_ARMCHAIR_BASE,
self.edgetype2, self.ncell1, growindex1)
# left and right lattice sides
elif growindex1 != self.ncell1 - 1 and growindex2 == 0 or \
growindex2 == self.ncell2 - 1:
atoms_to_terminate = get_edge_termination(
self.ZIGZAG_BY_ONE_BASE, self.ARMCHAIR_BY_ONE_BASE,
self.edgetype1, self.ncell2, growindex2)
return atoms_to_terminate
[docs]class HoneycombCell(object):
"""
Manage cells in the honeycomb lattice.
"""
[docs] def __init__(self, unit_cell_structure, center, growindex1, growindex2):
"""
Create an instance.
:type unit_cell_structure: schrodinger.structure.Structure
:param unit_cell_structure: unit cell structure
:type center: numpy.array
:param center: center of cell
:type growindex1: int
:param growindex1: index along first grow vector
:type growindex2: int
:param growindex2: index along second grow vector
"""
self.structure = unit_cell_structure.copy()
self.center = center
self.growindex1 = growindex1
self.growindex2 = growindex2
self.neighbors = {}
self.fragment = None
self.frag_to_latt = {}
[docs] def createCell(self):
"""
Make a cell.
"""
transform.translate_structure(self.structure, self.center[0],
self.center[1], self.center[2])
[docs] def findNeighboringCells(self, cells, deltas):
"""
For the current cell find the neighboring cells.
:type cells: OrderedDict
:type cells: keys are tuples of grow indicies and values are
HoneycombCell objects
:type deltas: list of tuples
:type deltas: contains deltas to potential neighboring cells
"""
# this function marks the difference between what is meant by
# potential neighbors, which are positions where a neighbor
# might be depending on the lattice shape, and neighbors,
# which are the actual neighbors, found by comparing with
# the current lattice structure in the middle of its growth
for delta in deltas:
indextuple = (self.growindex1 + delta[0],
self.growindex2 + delta[1])
neighboringcell = cells.get(indextuple)
if neighboringcell:
self.neighbors[delta] = neighboringcell
[docs] def makeFragment(self, redundantatoms):
"""
Create the cell fragment to be merged with the honeycomb lattice.
:type redundantatoms: list of ints
:param redundantatoms: indicies of redundant atoms present in the cell
"""
self.fragment = self.structure.copy()
self.fragment.deleteAtoms(redundantatoms)
[docs] def determineLatticeNumbering(self, latticelen):
"""
Determine the map of cell fragment unit cell atom indicies to lattice
atom indicies.
:type latticelen: int
:param latticelen: current number of atoms in lattice
"""
for index, atom in enumerate(self.fragment.atom, 1):
fragindex = atom.property[HoneycombUnitCell.UNITCELLINDEX]
latticeindex = latticelen + index
self.frag_to_latt[fragindex] = latticeindex
[docs]class HoneycombLattice(object):
r"""
Create the honeycomb lattice according to the following coordinate system::
----------------------------------------------------------------------
| Y |
| -------------------- |
| / \ |
| / \ |
| / \ |
| / \ |
| / \ |
| / \ |
| -------------------- P2 | |
| / 6 5 \ / |
| / \ / |
| / \ / |
| / \ / |
| / \ / |
| / \ / |
| | 1 O 4 -------------------- X |
| \ / \ |
| \ / \ |
| \ / \ |
| \ / \ |
| \ / \ |
| \ 2 3 / \ |
| -------------------- P1 | |
| \ / |
| \ / |
| \ / |
| \ / |
| \ / |
| \ / |
| -------------------- |
----------------------------------------------------------------------
where in addition to the parameters specified in the HoneycombUnitCell class
the two lattice vectors are given by vec(P1, O) and vec(P2, O) with
angle(P1, O, P2) = 60 degrees.
"""
CAPPINGELEMENT = 'H'
CG_CAPPING_ELEMENT = 'C'
CAPPING_BOND_LENGTH = 1.103
DIRECTION = 'forward'
TITLEKEY = 's_m_title'
ENTRYKEY = 's_m_entry_name'
TITLENAME = 'nanosheet'
NCELL1 = 'i_matsci_N_Cell_1'
NCELL2 = 'i_matsci_N_Cell_2'
C_VACUUM = constants.Constants.BILAYERSEP
[docs] def __init__(self,
element1,
element2,
bondlength,
ncell1,
edgetype1,
ncell2,
edgetype2,
termfrag,
min_term_frags,
is_coarse_grain=False):
"""
Create an instance.
:type element1: str
:param element1: elemental symbol of the first atom
:type element2: str
:param element2: elemental symbol of the second atom
:type bondlength: float
:param bondlength: bond length between the first and
second atoms in Angstrom
:type ncell1: int
:param ncell1: number of cells along lattice side 1
:type edgetype1: str
:param edgetype1: type of edge for lattice side 1
:type ncell2: int
:param ncell2: number of cells along lattice side 2
:type edgetype2: str
:param edgetype2: type of edge for lattice side 2
:type termfrag: str
:param termfrag: terminate the lattice with a given fragment
:type min_term_frags: bool
:param min_term_frags: minimize the geometry of terminating
fragments
:type is_coarse_grain: bool
:param is_coarse_grain: Whether a coarse grain structure is being
created
"""
self.is_coarse_grain = is_coarse_grain
self.element1 = element1
self.element2 = element2
self.bondlength = bondlength
self.ncell1 = ncell1
self.edgetype1 = edgetype1
self.ncell2 = ncell2
self.edgetype2 = edgetype2
self.termfrag = termfrag
self.min_term_frags = min_term_frags
self.cells = OrderedDict()
self.structure = None
self.terminatingatoms = []
self.frozenatoms = []
self.atomic_number1 = None
self.atomic_number2 = None
[docs] def updateLatticeBonding(self, cell, bonds_to_neighbors):
"""
Update the bonding between the lattice and the newly added fragment.
:type cell: HoneycombCell
:param cell: recently added cell
:type bonds_to_neighbors: dict
:param bonds_to_neighbors: keys are tuples to neighboring cells
and values are lists of bonding atom pairs using the unit cell
numbering
"""
# loop over neighbors, get neighbor-cell bonding atom pairs,
# get lattice numbering
for delta, neighbor in cell.neighbors.items():
bondstoform = bonds_to_neighbors.get(delta)
for bond in bondstoform:
neighborindex = neighbor.frag_to_latt.get(bond[0])
cellindex = cell.frag_to_latt.get(bond[1])
if neighborindex and cellindex:
self.structure.addBond(self.structure.atom[neighborindex],
self.structure.atom[cellindex], 1)
self.structure.retype()
[docs] def terminateLattice(self, logger=None):
"""
Terminate the lattice with the given fragments.
:type logger: logging.getLogger
:param logger: output logger
"""
# this is done explicity rather than by the structutils.build.add_hydrogens()
# method because by default amino hydrogens are not added in a planar
# fashion and by default boron will be sp3, i.e. two hydrogens will be
# added. While we are doing the termination extend the list of frozen atoms
# to include the first atom in the terminating fragment, i.e. the one
# bonded to the terminating atom
capping_element = (self.CG_CAPPING_ELEMENT
if self.is_coarse_grain else self.CAPPINGELEMENT)
capping_bond_length = (self.bondlength if self.is_coarse_grain else
self.CAPPING_BOND_LENGTH)
for index in self.terminatingatoms:
terminatingatom = self.structure.atom[index]
# form vectors to atoms bonded to the terminating atom
atom1, atom2 = terminatingatom.bonded_atoms
bonded_atoms = set([atom1.index, atom2.index])
atom1vec = numpy.array(atom1.xyz) - numpy.array(terminatingatom.xyz)
atom2vec = numpy.array(atom2.xyz) - numpy.array(terminatingatom.xyz)
# the sum of these two vectors scaled provides the vector where the first
# fragment atom would be placed
new_atom_vec = transform.get_normalized_vector(atom1vec + atom2vec)
new_atom_vec = -1 * capping_bond_length * new_atom_vec
# get new vectors absolute position and add and bond the hydrogen
# (to be replaced by a non-hydrogen fragment if requested)
new_atom_vec = numpy.array(terminatingatom.xyz) + new_atom_vec
hatom = self.structure.addAtom(capping_element, new_atom_vec[0],
new_atom_vec[1], new_atom_vec[2])
if self.is_coarse_grain:
# Assume the covalent radius is approximately equal to the
# radius of the coarse grain bead
set_atom_coarse_grain_properties(self.structure,
hatom,
self.termfrag,
radius=capping_bond_length / 2)
frozenatom = hatom.index
self.structure.addBond(terminatingatom, hatom, 1)
# attach non-hydrogen fragment
if (self.termfrag != constants.Constants.TERMFRAG and
not self.is_coarse_grain):
fraggroup, fragname = constants.Constants.TERMDICT[
self.termfrag]
build.attach_fragment(self.structure, terminatingatom.index,
hatom.index, fraggroup, fragname,
self.DIRECTION)
new_terminatingatom = self.structure.atom[index]
new_bonded_atoms = \
set([aatom.index for aatom in new_terminatingatom.bonded_atoms])
frozenatom = new_bonded_atoms.difference(bonded_atoms).pop()
self.frozenatoms.append(frozenatom)
self.structure.retype()
color.apply_color_scheme(self.structure, 'element')
[docs] def minTerminatingFragments(self, logger=None):
"""
Minimize the geometry of all atoms in the terminating fragments with
the exception of the lattice bound atoms.
:type logger: logging.getLogger
:param logger: output logger
"""
# minimize everything other than the lattice itself as well
# as those atoms bound directly to the lattice. The minimization
# might really only need to be done over the rotatable bond
# connecting fragment to lattice, however, that appears to
# currently be more trouble than it is worth considering
# that the minimize.Minimizer class only allows restraints.
with minimize.minimizer_context() as minimizeobj:
try:
minimizeobj.setStructure(self.structure)
except ValueError:
if logger:
msg = ('The structure can not be minimized. The '
'unminimized structure will be provided.')
logger.warning(msg)
return
for frozenatom in self.frozenatoms:
minimizeobj.addPosFrozen(frozenatom)
minimizeobj.minimize()
[docs] def setLatticeProperties(self, chorus_properties):
"""
Set some structure properties of the lattice.
:type chorus_properties: list
:param chorus_properties: contains the nine chorus properties,
i.e. ax, ay, az, bx, ..., cz
"""
for atom in self.structure.atom:
try:
del atom.property[HoneycombUnitCell.UNITCELLINDEX]
except KeyError:
# It is ok if some atoms do not have this property
pass
self.structure.property[self.TITLEKEY] = self.TITLENAME
self.structure.property[self.ENTRYKEY] = self.TITLENAME
self.structure.property[self.NCELL1] = self.ncell1
self.structure.property[self.NCELL2] = self.ncell2
# if not terminating then set up PBC
if self.termfrag == constants.Constants.TERMFRAGS[0]:
xtal.set_pbc_properties(self.structure, chorus_properties)
self.structure.property[xtal.Crystal.SPACE_GROUP_KEY] = \
xtal.P1_SPACE_GROUP_SYMBOL
[docs] def getChorusPBC(self, growvec1, growvec2):
"""
Return the chorus box PBC.
:type growvec1: numpy.array
:param growvec1: the first lattice grow vector
:type growvec2: numpy.array
:param growvec2: the second lattice grow vector
:rtype: list
:return: contains the nine chorus properties, i.e. ax, ay, az, bx, ..., cz
"""
a_vec = growvec1 * self.ncell1
b_vec = growvec2 * self.ncell2
c_vec = numpy.array([0.0, 0.0, self.C_VACUUM])
return list(a_vec) + list(b_vec) + list(c_vec)
[docs] def buildLattice(self, logger=None):
"""
Build a honeycomb lattice.
:type logger: logging.getLogger
:param logger: output logger
"""
# this is the main function that orchestrates everything
# get a unit cell
unitcell = HoneycombUnitCell(self.element1, self.element2,
self.bondlength, self.is_coarse_grain)
unitcell.buildUnitCell()
self.atomic_number1 = unitcell.atomic_number1
self.atomic_number2 = unitcell.atomic_number2
# determine the lattice vectors
lattvec1, lattvec2 = get_lattice_vectors(self.bondlength)
# determine grow parameters, i.e. grow vectors, connections to
# neighbors, etc.
grow = Grow(self.edgetype1, self.edgetype2, self.ncell1, self.ncell2)
growvec1, growvec2 = grow.getGrowVectors(lattvec1, lattvec2)
neighbortypes = grow.getNeighborTypes()
# make lattice row-wise
for growindex1 in range(self.ncell1):
for growindex2 in range(self.ncell2):
# make a cell with a certain center, keep a list
# of added cells
center = growindex1 * growvec1 + growindex2 * growvec2
cell = HoneycombCell(unitcell.structure, center, growindex1,
growindex2)
self.cells[(growindex1, growindex2)] = cell
cell.createCell()
# handle first cell in lattice as a special case, for
# other cells (1) find neighbors, (2) remove redundant
# atoms in current cell, (3) add fragment to lattice,
# (4) bond fragment to lattice
if not growindex1 and not growindex2:
cell.fragment = cell.structure.copy()
indicies = list(range(1, HoneycombUnitCell.NUMATOMS + 1))
cell.frag_to_latt = dict(list(zip(indicies, indicies)))
self.structure = cell.fragment.copy()
else:
cell.findNeighboringCells(self.cells, list(neighbortypes))
redundantatoms = grow.getRedundantAtoms(
list(cell.neighbors), neighbortypes)
cell.makeFragment(redundantatoms)
cell.determineLatticeNumbering(len(self.structure.atom))
self.structure.extend(cell.fragment)
bonds_to_neighbors = \
grow.getBondingInfo(list(cell.neighbors), neighbortypes)
self.updateLatticeBonding(cell, bonds_to_neighbors)
# update list of lattice atoms to terminate
atoms_to_terminate = \
grow.getAtomsToTerminate(growindex1, growindex2, logger)
terminatingatoms = [cell.frag_to_latt.get(atom) \
for atom in atoms_to_terminate]
self.terminatingatoms.extend(terminatingatoms)
# terminate the lattice
if self.termfrag != constants.Constants.TERMFRAGS[0]:
self.frozenatoms.extend(list(range(1,
len(self.structure.atom) + 1)))
self.terminateLattice(logger)
if len(self.frozenatoms) != len(self.structure.atom) and \
self.min_term_frags:
self.minTerminatingFragments(logger)
chorus_properties = self.getChorusPBC(growvec1, growvec2)
self.setLatticeProperties(chorus_properties)
[docs]def get_lattice_vectors(bondlength):
"""
Return the two honeycomb lattice vectors for a nanosheet. We build the
lattice in XY-plane, where Z = 0
:param float bondlength: The bond length between the atoms comprising the
nanosheet (an Angstroms)
:return (numpy.array, numpy.array): A 2-tuple of the two lattice vectors
"""
length = 2 * bondlength * math.sin(HoneycombUnitCell.ANGLEMEDIUM)
x = length * math.cos(HoneycombUnitCell.ANGLESMALL)
y = length * math.sin(HoneycombUnitCell.ANGLESMALL)
z = float(0)
lattvec1 = numpy.array([x, -1 * y, z])
lattvec2 = numpy.array([x, y, z])
return lattvec1, lattvec2
[docs]class HoneycombBilayers(object):
"""
Create honeycomb bilayers.
"""
ZEROVEC = numpy.zeros(3)
STACK_DIRECTION = -1
[docs] def __init__(self, lattice, separation, bondlength, stacktype, nbilayers,
bilayershift):
"""
Create an instance.
:type lattice: schrodinger.structure.Structure
:param lattice: lattice structure
:type separation: float
:param separation: bilayer separation
:type bondlength: float
:param bondlength: bond length between the first and
second atoms in Angstrom
:type stacktype: str
:param stacktype: type of bilayer stacking to be used
:type nbilayers: int
:param nbilayers: number of bilayers
:type bilayershift: float
:param bilayershift: offset of bilayers in terms of
the number of unit cells
"""
self.lattice = lattice
self.separation = separation
self.bondlength = bondlength
self.stacktype = stacktype
self.nbilayers = nbilayers
self.bilayershift = bilayershift
# these are the bilayer stacking vectors
self.xvec = 2 * self.bilayershift * self.bondlength * numpy.array(
transform.X_AXIS)
self.zvec = self.STACK_DIRECTION * self.separation * numpy.array(
transform.Z_AXIS)
self.structure = lattice.copy()
self.buildBilayer()
self.stackBilayers()
[docs] def buildBilayer(self):
"""
Build the bilayer.
"""
# copy the lattice and translate it from the original
xvec = self.xvec
zvec = self.zvec
transvec = xvec + zvec
lattice = self.lattice.copy()
transform.translate_structure(lattice, transvec[0], transvec[1],
transvec[2])
self.structure.extend(lattice)
[docs] def stackBilayers(self):
"""
Stack the bilayers.
"""
# copy the bilayer and translate it from the original
bilayer = self.structure.copy()
xvec = self.ZEROVEC
for bilayerindex in range(2, self.nbilayers + 1):
if self.stacktype == constants.Constants.ABCD:
xvec = 2 * (bilayerindex - 1) * self.xvec
zvec = 2 * (bilayerindex - 1) * self.zvec
transvec = xvec + zvec
nextbilayer = bilayer.copy()
transform.translate_structure(nextbilayer, transvec[0], transvec[1],
transvec[2])
self.structure.extend(nextbilayer)
[docs] def updatePBC(self):
"""
Update the PBC.
"""
a_vec, b_vec, c_vec = xtal.get_vectors_from_chorus(self.structure)
c_vec = 2 * self.STACK_DIRECTION * self.nbilayers * self.zvec
chorus_properties = list(a_vec) + list(b_vec) + list(c_vec)
xtal.set_pbc_properties(self.structure, chorus_properties)
[docs] def translateLayers(self):
"""
Translate the layers to be inside the box.
"""
transvec = (2 * self.nbilayers - 1) * self.STACK_DIRECTION * self.zvec
transform.translate_structure(self.structure, *transvec)
[docs] def getTerminatingAtoms(self, terminatingatoms):
"""
Return the terminating atoms for all layers.
:type terminatingatoms: list
:param termiatingatoms: terminating atoms for the first layer
:rtype: list
:return: terminating atoms for all layers
"""
natoms = self.lattice.atom_total
all_terminatingatoms = list(terminatingatoms)
for i_idx in range(2, 2 * self.nbilayers + 1):
indices = [
j_idx + (i_idx - 1) * natoms for j_idx in terminatingatoms
]
all_terminatingatoms.extend(indices)
return all_terminatingatoms
[docs]class NanoSheets(object):
"""
Main class for making nanosheets.
"""
MSGWIDTH = 50
[docs] def __init__(self,
element1=constants.Constants.ELEMENT1,
element2=constants.Constants.ELEMENT2,
bondlength=constants.Constants.BONDLENGTH,
ncell1=constants.Constants.NCELL1,
edgetype1=constants.Constants.EDGETYPE1,
ncell2=constants.Constants.NCELL2,
edgetype2=constants.Constants.EDGETYPE2,
no_double_bonds=constants.Constants.NO_DOUBLE_BONDS,
termfrag=constants.Constants.TERMFRAG,
min_term_frags=constants.Constants.MIN_TERM_FRAGS,
bilayersep=constants.Constants.BILAYERSEP,
nbilayers=constants.Constants.NBILAYERS,
stacktype=constants.Constants.ABAB,
bilayershift=constants.Constants.BILAYERSHIFT,
orient=False,
logger=None,
is_coarse_grain=False):
"""
:type element1: str
:param element1: elemental symbol of the first atom
:type element2: str
:param element2: elemental symbol of the second atom
:type bondlength: float
:param bondlength: bond length between the first and
second atoms in Angstrom
:type ncell1: int
:param ncell1: number of cells along lattice side 1
:type edgetype1: str
:param edgetype1: type of edge for lattice side 1
:type ncell2: int
:param ncell2: number of cells along lattice side 2
:type edgetype2: str
:param edgetype2: type of edge for lattice side 2
:type no_double_bonds: bool
:param no_double_bonds: disable the formation of double bonds
:type termfrag: str
:param termfrag: terminate the lattice with a given fragment
:type min_term_frags: bool
:param min_term_frags: minimize the geometry of terminating
fragments
:type bilayersep: float
:param bilayersep: bilayer separation
:type nbilayers: int
:param nbilayers: number of bilayers
:type stacktype: str
:param stacktype: bilayer stacking type
:type bilayershift: float
:param bilayershift: offset of bilayer in terms of the
number of unit cells
:type orient: bool
:param orient: whether to orient the sheets for Maestro
:type logger: logging.getLogger
:param logger: output logger
:type is_coarse_grain: bool
:param is_coarse_grain: Whether a coarse grain structure is being
created
"""
# check user options
chkobj = CheckInput()
chkobj.checkAll(element1, element2, bondlength, ncell1, edgetype1,
ncell2, edgetype2, no_double_bonds, termfrag,
min_term_frags, bilayersep, nbilayers, stacktype,
bilayershift, logger, is_coarse_grain)
self.element1 = chkobj.element1
self.element2 = chkobj.element2
self.bondlength = chkobj.bondlength
self.ncell1 = chkobj.ncell1
self.edgetype1 = chkobj.edgetype1
self.ncell2 = chkobj.ncell2
self.edgetype2 = chkobj.edgetype2
self.no_double_bonds = chkobj.no_double_bonds
self.termfrag = chkobj.termfrag
self.min_term_frags = chkobj.min_term_frags
self.bilayersep = chkobj.bilayersep
self.nbilayers = chkobj.nbilayers
self.stacktype = chkobj.stacktype
self.bilayershift = chkobj.bilayershift
self.orient = orient
self.is_coarse_grain = is_coarse_grain
self.structure = None
# print job parameters
if logger:
self.printJobParams(logger)
# make sheets
self.makeNanoSheets(logger)
[docs] def printJobParams(self, logger=None):
"""
Print job parameters.
:type logger: logging.getLogger
:param logger: output logger
"""
# build formatted strings first
element1 = tlog.get_param_string('Element 1', self.element1,
self.MSGWIDTH)
element2 = tlog.get_param_string('Element 2', self.element2,
self.MSGWIDTH)
bondlength = tlog.get_param_string('Bond length / Ang.',
self.bondlength, self.MSGWIDTH)
no_double_bonds = tlog.get_param_string('No double bonds',
self.no_double_bonds,
self.MSGWIDTH)
ncell1 = tlog.get_param_string('Number of cells along side 1',
self.ncell1, self.MSGWIDTH)
ncell2 = tlog.get_param_string('Number of cells along side 2',
self.ncell2, self.MSGWIDTH)
edgetype1 = tlog.get_param_string('Edge type along side 1',
self.edgetype1, self.MSGWIDTH)
edgetype2 = tlog.get_param_string('Edge type along side 2',
self.edgetype2, self.MSGWIDTH)
termfrag = tlog.get_param_string('Terminating fragment', self.termfrag,
self.MSGWIDTH)
min_term_frags = tlog.get_param_string('Minimize terminating fragments',
self.min_term_frags,
self.MSGWIDTH)
nbilayers = tlog.get_param_string('Number of bilayers', self.nbilayers,
self.MSGWIDTH)
bilayersep = tlog.get_param_string('Bilayer separation / Ang.',
self.bilayersep, self.MSGWIDTH)
bilayershift = tlog.get_param_string('Bilayer shift', self.bilayershift,
self.MSGWIDTH)
stacktype = tlog.get_param_string('Bilayer stacking type',
self.stacktype, self.MSGWIDTH)
SEPARATOR = self.MSGWIDTH * '='
# print formatted strings
logger.info('Job Parameters:')
logger.info('')
logger.info(SEPARATOR)
logger.info('')
logger.info('Unit cell')
logger.info('---------')
logger.info('')
logger.info(element1)
logger.info(element2)
logger.info(bondlength)
logger.info(no_double_bonds)
logger.info('')
logger.info('Lattice')
logger.info('-------')
logger.info('')
logger.info(ncell1)
logger.info(ncell2)
logger.info(edgetype1)
logger.info(edgetype2)
logger.info(termfrag)
logger.info(min_term_frags)
logger.info('')
if self.nbilayers:
logger.info('Bilayer')
logger.info('-------')
logger.info('')
logger.info(nbilayers)
logger.info(bilayersep)
logger.info(bilayershift)
logger.info(stacktype)
logger.info('')
logger.info(SEPARATOR)
logger.info('')
[docs] def makeNanoSheets(self, logger=None):
"""
Make nanosheets.
:type logger: logging.getLogger
:param logger: output logger
"""
lattice = HoneycombLattice(self.element1, self.element2,
self.bondlength, self.ncell1, self.edgetype1,
self.ncell2, self.edgetype2, self.termfrag,
self.min_term_frags, self.is_coarse_grain)
lattice.buildLattice(logger)
self.structure = lattice.structure
terminatingatoms = lattice.terminatingatoms
if self.nbilayers:
bilayers = HoneycombBilayers(lattice.structure, self.bilayersep,
self.bondlength, self.stacktype,
self.nbilayers, self.bilayershift)
if self.termfrag == constants.Constants.TERMFRAGS[0]:
bilayers.updatePBC()
if self.orient:
bilayers.translateLayers()
terminatingatoms = bilayers.getTerminatingAtoms(
lattice.terminatingatoms)
self.structure = bilayers.structure
if self.orient and self.termfrag == constants.Constants.TERMFRAGS[0]:
vectors = xtal.get_vectors_from_chorus(self.structure)
self.structure, origin, a_vec, b_vec, c_vec = \
slab.maestro_rotate_cell(self.structure, numpy.array(xtal.ParserWrapper.ORIGIN), \
*vectors)
chorus_properties = list(a_vec) + list(b_vec) + list(c_vec)
xtal.set_pbc_properties(self.structure, chorus_properties)
self.structure = build_cell(self.structure, terminatingatoms,
lattice.atomic_number1,
lattice.atomic_number2, self.bondlength)
if self.is_coarse_grain:
# Note do not set this property during build as it changes the
# covalent radius calculations
set_as_coarse_grain(self.structure)
elif not self.no_double_bonds:
self.structure = xtal.assign_bond_orders_w_mmlewis(self.structure,
fix_metals=False,
logger=logger)
[docs]def build_cell(astructure,
terminatingatoms,
atomic_number1,
atomic_number2,
bondlength,
no_bonding_along=None):
"""
Build the cell.
:type astructure: schrodinger.structure.Structure
:param astructure: the structure for which to build the cell
:type terminatingatoms: list
:param terminatingatoms: the terminating atoms for which PBC bonds
will be needed
:type atomic_number1: int
:param atomic_number1: the atomic number of the first element
:type atomic_number2: int
:param atomic_number2: the atomic number of the second element
:type bondlength: float
:param bondlength: the bond length
:type no_bonding_along: list or None
:param no_bonding_along: contains the indices of the lattice vectors along
which to not form bonds, indices are 0 (a-vector), 1 (b-vector), and 2
(c-vector), if None then [2] will be used which is the default for sheets
:rtype: schrodinger.structure.Structure
:return: the built cell
"""
PRIME_LEN = 1000.
BOND_THRESH = 0.01
# delete duplicate atoms and determine new terminating atoms
tmp_st = astructure.copy()
renumber_map = xtal.delete_duplicate_atoms(tmp_st)
if not renumber_map:
renumber_map = {idx: idx for idx in range(1, tmp_st.atom_total + 1)}
new_terminatingatoms = []
for oidx, nidx in renumber_map.items():
if nidx is None:
aatom = astructure.atom[oidx]
for batom in aatom.bonded_atoms:
if renumber_map[batom.index]:
new_terminatingatoms.append(batom.index)
astructure = tmp_st
# update terminating atoms
terminatingatoms = [renumber_map[idx] for idx in terminatingatoms \
if renumber_map[idx]]
terminatingatoms.extend([renumber_map[idx] for idx in new_terminatingatoms \
if renumber_map[idx]])
terminatingatoms.sort()
# this temporarily increases the lengths of the lattice vectors
# along which there will be no bonds formed, i.e. for sheets the
# c-vector and for tubes the a- and c-vectors, this is to allow
# connect_atoms to form PBC bonds for cases where the user
# has specified large bond lengths
if no_bonding_along is None:
no_bonding_along = [2]
orig_vecs = numpy.array(xtal.get_vectors_from_chorus(astructure))
new_vecs = numpy.copy(orig_vecs)
for idx in no_bonding_along:
new_vecs[idx] = PRIME_LEN * transform.get_normalized_vector(
new_vecs[idx])
chorus_properties = list(new_vecs.flatten())
xtal.set_pbc_properties(astructure, chorus_properties)
# shift the covalent bonding parameters to honor the requested
# bond length and define the PBC bonds
cov_min = bondlength - BOND_THRESH
cov_radii = xtal.get_cov_radii()
cov_rad_a = cov_radii[atomic_number1]
cov_rad_b = cov_radii[atomic_number2]
cov_offset = bondlength + BOND_THRESH - cov_rad_a - cov_rad_b
cov_factor = xtal.ParserWrapper.COV_FACTOR
xtal.connect_atoms(astructure,
atoms_to_connect=terminatingatoms,
cov_min=cov_min,
cov_offset=cov_offset,
cov_factor=cov_factor,
delete_existing=False,
pbc_bonding=True,
only_pbc_bonds=True)
# reinstate the original c vector length
chorus_properties = list(orig_vecs.flatten())
xtal.set_pbc_properties(astructure, chorus_properties)
# build the cell using the existing PBC bonds
CHOICE_OFF = xtal.ParserWrapper.CHOICE_OFF
XTAL_DICT = dict(
list(zip(['bonding', 'bond_orders', 'translate'], 3 * [CHOICE_OFF])))
xtal_kwargs = dict(XTAL_DICT)
xtal_kwargs['use_existing_pbc_bonds'] = True
astructure = xtal.get_cell(astructure, xtal_kwargs=xtal_kwargs)
# MATSCI-3173 - by default translate atoms to first unit cell
xtal.translate_to_cell(astructure)
xtal.label_pbc_bonds(astructure)
return astructure