"""
Classes and functions for creating crystals by unit cell.
Copyright Schrodinger, LLC. All rights reserved."""
# Contributor: Thomas F. Hughes
import argparse
from collections import namedtuple
import copy
import itertools
import math
import os
import re
import warnings
from functools import reduce
from past.utils import old_div
import numpy
import spglib
from scipy import constants
import pymmlibs
import schrodinger.structure as structure
from schrodinger.application.matsci import coarsegrain
from schrodinger.application.matsci import msprops
from schrodinger.application.matsci import textlogger
from schrodinger.application.matsci.nano import space_groups
from schrodinger.infra import mm
from schrodinger.infra import structure as infrastructure
from schrodinger.infra.mm import \
mmelement_get_atomic_number_by_symbol as get_atomic_number
from schrodinger.infra.mm import \
mmelement_get_atomic_weight_by_atomic_number as get_atomic_weight
from schrodinger.infra.mmerr import ErrorHandler
from schrodinger.structutils import analyze
from schrodinger.structutils import assignbondorders
from schrodinger.structutils import build
from schrodinger.structutils import transform
from schrodinger.test import ioredirect
from schrodinger.utils import fileutils
from schrodinger.utils import mmutil
from schrodinger.utils import subprocess
from schrodinger.application.desmond.constants import CT_TYPE
from schrodinger.application.matsci import clusterstruct
from schrodinger.application.matsci import desmondutils # noqa: F401
from schrodinger.application.matsci import elementalprops
from schrodinger.application.matsci import msutils
_version = '$Revision 0.0 $'
gcd = math.gcd
# Space group property name
SPACE_GROUP_KEY = mm.M2IO_PDB_CRYSTAL_SPACE_GROUP
# Space group ID
SPACE_GROUP_ID_KEY = msprops.SPACE_GROUP_ID_KEY
# P1 space group symbol
P1_SPACE_GROUP_SYMBOL = mm.P1_SPACE_GROUP
P1_SPACE_GROUP_ID = 1
COV_RADIUS_KEY = 'r_matsci_COV_radius/Ang.'
# this threshold is used to pin floating point fractional
# coordinates and in particular the ones on cell boundaries
FRACT_OFFSET = 0.0001
# Interatomic distance for which atoms are considered to occupy the same
# position (overlapping)
OVERLAP_ATOM_THRESHOLD = 0.25
# the following list is needed to determine the possible
# tail-to-head directions of PBC bonds in cells in a super
# cell, meaning given a PBC bond in a cell in a super cell,
# which cells of the super cell can that PBC bond point
# ahead and behind to, of the twenty-six, six along faces,
# twelve along edges, and eight along corners, spanning
# vectors for a cell, we only need half of them, it is
# easier to enumerate them rather than use
# itertools.product([-1, 0, 1], repeat=3) and then trim,
# because we grow the super cell along the c direction 1st,
# then along the b direction 2nd, and then along the a
# direction 3rd, it is important to push the negative
# coefficients to c firstly and b secondly this way moving
# behind a cell will point firstly along c and secondly
# along b, and ahead in the opposite way
AHEAD_BASE_TRIPLES = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
AHEAD_TRIPLES = AHEAD_BASE_TRIPLES + [(0, 1, 1), (0, 1, -1), (1, 0, 1),
(1, 1, 0), (1, 1, 1), (1, 1, -1),
(1, 0, -1), (1, -1, 0), (1, -1, -1),
(1, -1, 1)]
MAX_BONDS_EXCEPTIONS = {'Pr': 10, 'O': 5, 'Na': 6, 'Cl': 6, 'I': 6}
LATTICE_PARAMS_KEYS = ['a_param', 'b_param', 'c_param', \
'alpha_param', 'beta_param', 'gamma_param']
EXTENTS_KEYS = ['ncella', 'ncellb', 'ncellc']
# checked these combinations by hand to make sure they make sense
# all look normal
ATOM_TO_BOND_REPR_DICT = {
structure.ATOM_NOSTYLE: structure.BOND_WIRE,
structure.ATOM_CIRCLE: structure.BOND_WIRE,
structure.ATOM_CPK: structure.BOND_TUBE,
structure.ATOM_BALLNSTICK: structure.BOND_BALLNSTICK
}
PBC_POSITION_KEY = msprops.PBC_POSITION_KEY
CENTER_PBC_POSITION = 'center_structure'
ANCHOR_PREFIX = 'anchor_'
# slots are x, y, and z in Cartesian coordinates
ANCHOR_PBC_POSITION = ANCHOR_PREFIX + '%s_%s_%s'
# Doubly named space groups in spglib
SPGLIB_DOUBLY_NAMED_SHORT = {'Pncb', 'Pcna', 'Pnmm', 'Pmnm', 'Bbeb', 'Bbcb'}
# Default precision for space group assignment
ASSIGN_SPG_SYMPREC = 1e-2
# Precision thresholds used in get_normal_surf_from_HKL
NORMAL_SURF_HKL_THR = 1e-8
SPG_TO_HALL_NUMBER = [
1, 2, 3, 6, 9, 18, 21, 30, 39, 57, 60, 63, 72, 81, 90, 108, 109, 112, 115,
116, 119, 122, 123, 124, 125, 128, 134, 137, 143, 149, 155, 161, 164, 170,
173, 176, 182, 185, 191, 197, 203, 209, 212, 215, 218, 221, 227, 228, 230,
233, 239, 245, 251, 257, 263, 266, 269, 275, 278, 284, 290, 292, 298, 304,
310, 313, 316, 322, 334, 335, 337, 338, 341, 343, 349, 350, 351, 352, 353,
354, 355, 356, 357, 358, 359, 361, 363, 364, 366, 367, 368, 369, 370, 371,
372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386,
387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401,
402, 404, 406, 407, 408, 410, 412, 413, 414, 416, 418, 419, 420, 422, 424,
425, 426, 428, 430, 431, 432, 433, 435, 436, 438, 439, 440, 441, 442, 443,
444, 446, 447, 448, 449, 450, 452, 454, 455, 456, 457, 458, 460, 462, 463,
464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478,
479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493,
494, 495, 497, 498, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510,
511, 512, 513, 514, 515, 516, 517, 518, 520, 521, 523, 524, 525, 527, 529,
530] # yapf: disable
LATT_ROUND = 6
LAST_ELEMENTS = ('H', 'C', 'N', 'O', 'F', 'P', 'S', 'Cl', 'Se', 'Br', 'I', 'At')
_COV_RADII = None
# Maximum cell length to enable check check_also_reg_bond in is_pbc_bond
SMALL_CELL_LENGTH = 10.
PBC_BUFFER_LENGTH = 2 # unit: (A)
AXIS = namedtuple('AXIS', 'X Y Z', defaults=(0, 1, 2))
[docs]def find_furthest_atom_along_vector(struct, vector):
"""
Find the atom in struct that is furthest in the direction of vector. For
instance, if vector is the Z-axis, this will find the atom with the largest
z coordinate
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to check
:type vector: `numpy.array`
:param vector: The vector that gives the direction to search
:rtype: `schrodinger.structure._StructureAtom`
:return: The atom furthest in the vector direction
"""
maxdist = -numpy.inf
dist_atom = None
for atom in struct.atom:
# The atom whose coordinate have the largest dot product with vector is
# the furthest atom
distance = vector.dot(atom.xyz)
if distance > maxdist:
maxdist = distance
dist_atom = atom
return dist_atom
[docs]def find_origin_on_structure_exterior(struct, vector):
"""
Find the point to originate the normal plane to the passed vector so that
the passed structure is entirely behind the vector (i.e. the structure will
be entirely on one side of the plane, and vector on the other side).
"""
# Find the atom with the largest component in the direction of the
# plane normal. This is the atom that defines the location of the
# plane along the plane normal direction.
dist_atom = find_furthest_atom_along_vector(struct, vector)
vorigin = numpy.array(dist_atom.xyz)
# The vector origin is now at the atom furthest in the vector's
# direction. Move the origin of the vector so that it is directly over the
# center of mass of the structure
vector_to_centroid = transform.get_centroid(struct)[:3] - vorigin
# Remove the component of the atom->origin vector that is parallel to
# the plane normal - that will leave us a vector that slides along the
# plane.
normvec = transform.get_normalized_vector(vector)
slide = vector_to_centroid - vector_to_centroid.dot(normvec) * normvec
# Now slide the vector along the plane
vorigin = vorigin + slide
return vorigin
[docs]def is_infinite(astructure):
"""
Return a boolean indicating whether the given structure
is infinitely bonded, meaning that it has bonds that
cross the periodic boundary which can not be eliminated
by means of molecule contraction. If any molecule in the
given structure is infinitely bonded then the structure
itself is considered infinite. As examples of infinite
systems consider graphene, gold, infinite polymer, etc.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure whose infiniteness is
in question, can be a unit cell or an ASU
:raise ValueError: if there is an issue with the input
:rtype: bool
:return: True if infinite, False otherwise
"""
try:
lattice_params = get_lattice_param_properties(astructure)
except KeyError:
lattice_params = None
space_group = astructure.property.get(Crystal.SPACE_GROUP_KEY)
if lattice_params and not space_group or not lattice_params and space_group:
msg = ('Set of lattice PDB properties, i.e. '
'*_pdb_PDB_CRYST1_*, are incomplete.')
raise ValueError(msg)
try:
chorus_properties = get_chorus_properties(astructure)
except KeyError:
chorus_properties = None
if not chorus_properties and not lattice_params:
return False
build = not chorus_properties
if build:
on = ParserWrapper.CHOICE_ON
off = ParserWrapper.CHOICE_OFF
xtal_kwargs = {'bonding': on, 'bond_orders': off, 'translate': off}
cell = get_cell(astructure, xtal_kwargs=xtal_kwargs)
else:
cell = astructure.copy()
cell = sync_pbc(cell)
info = get_unit_lattice_vector_info(cell)
pbc = infrastructure.PBC(*get_chorus_properties(cell))
if not build:
for abond in cell.bond:
pbc_bond, also_reg_bond, dist = is_pbc_bond(
cell,
abond.atom1,
abond.atom2,
check_also_reg_bond=True,
unit_lattice_vectors=info,
pbc=pbc)
if pbc_bond:
if also_reg_bond:
return True
else:
clusterstruct.contract_structure(cell)
build = True
break
if build:
for abond in cell.bond:
pbc_bond, also_reg_bond, dist = is_pbc_bond(
cell,
abond.atom1,
abond.atom2,
check_also_reg_bond=True,
unit_lattice_vectors=info,
pbc=pbc)
if pbc_bond:
return True
return False
[docs]def get_check_also_reg_bond(struct, pbc=None, is_cg=None):
"""
Get the value of the check_also_reg_bond that is passed to is_pbc_bond based
on the PBC cell lengths.
:type astructure: `schrodinger.structure.Structure`
:param astructure: Structure for which to check the value of
check_also_reg_bond flag
:type pbc: `schrodinger.infra.structure.PBC` or None
:param pbc: The infrastructure PBC created from the Chorus box properties
:type is_cg: bool or None
:param is_cg: True, if structure is CG, otherwise False. If None, this will
be evaluated
:rtype: bool
:return: Flag value
"""
if is_cg is None:
is_cg = msutils.is_coarse_grain(struct)
if is_cg:
# get maximum particle radius in the CG structure
max_radius = max(
get_cov_radius(struct, atom.index, is_cg) for atom in struct.atom)
small_cell_length = 2. * max_radius
else:
small_cell_length = SMALL_CELL_LENGTH
if pbc is None:
pbc = infrastructure.PBC(*get_chorus_properties(struct))
return min(pbc.getBoxLengths()) < small_cell_length
[docs]def is_infinite2(astructure,
check_also_reg_bond=None,
is_cg=None,
atom_indices=None):
"""
Return a boolean indicating whether the given structure is infinitely
bonded, meaning that it has bonds that cross the periodic boundary which
can not be eliminated by means of molecule contraction. If any molecule in
the given structure is infinitely bonded then the structure itself is
considered infinite. As examples of infinite systems consider graphene,
gold, infinite polymer, etc. Note that, correct bond information is
expected in the input structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure whose infiniteness is
in question, can be a unit cell or an ASU
:type check_also_reg_bond: bool or None
:param check_also_reg_bond: check if the PBC bond is also a regular bond,
meaning that one of the atoms of the PBC bond is covalently bound to
two copies of the other atom, one inside the cell and one outside the
cell. If None, bool value will be decided automatically (preferred)
:type is_cg: bool
:param is_cg: True, if structure is CG, otherwise False. Needed if
check_also_reg_bond is None. If None, this will be evaluated
:type atom_indices: list or None
:paran atom_indices: Atom indices to check from the structure. If None,
entire stucture is checked
:raise KeyError: If chorus properties are missing
:rtype: bool
:return: True if infinite, False otherwise
"""
if atom_indices is None:
struct = astructure.copy()
else:
struct = astructure.extract(atom_indices, copy_props=True)
try:
chorus_properties = get_chorus_properties(struct)
except KeyError:
# Non-periodic struct is finite
return False
pbc = infrastructure.PBC(*chorus_properties)
info = get_unit_lattice_vector_info(struct)
if check_also_reg_bond is None:
if is_cg is None:
is_cg = msutils.is_coarse_grain(struct)
check_also_reg_bond = get_check_also_reg_bond(struct,
pbc=pbc,
is_cg=is_cg)
clusterstruct.contract_by_molecule(struct)
for mol in struct.molecule:
# This loop is similar to the one created in
# _StructureBondContainer::__iter__
for atom in mol.atom:
for abond in atom.bond:
if abond.atom2.index < atom.index:
continue
pbc_bond, also_reg_bond, dist = is_pbc_bond(
struct,
abond.atom1,
abond.atom2,
check_also_reg_bond=check_also_reg_bond,
unit_lattice_vectors=info,
pbc=pbc)
if pbc_bond:
return True
return False
[docs]def store_chorus_box_props(struct,
ax,
ay=0.0,
az=0.0,
bx=0.0,
by=None,
bz=0.0,
cx=0.0,
cy=0.0,
cz=None):
"""
Add properties to the structure that define the periodic boundary
condition in the way Desmond wants it defined.
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to add the properties to
:type ax: float
:param ax: The value of the ax box property
:type ay: float
:param ay: The value of the ay box property. Defaults to 0.
:type az: float
:param az: The value of the az box property. Defaults to 0.
:type bx: float
:param bx: The value of the bx box property. Defaults to 0.
:type by: float
:param by: The value of the by box property. If not given, this value is set
the same as ax.
:type bz: float
:param bz: The value of the bz box property. Defaults to 0.
:type cx: float
:param cx: The value of the cx box property. Defaults to 0.
:type cy: float
:param cy: The value of the cy box property. Defaults to 0.
:type cz: float
:param cz: The value of the cz box property. If not given, this value is set
the same as ax.
"""
if by is None:
by = ax
if cz is None:
cz = ax
# Float on the lines below is useful because sometimes these properties come
# from numpy arrays and are numpy double
struct.property[Crystal.CHORUS_BOX_AX_KEY] = float(ax)
struct.property[Crystal.CHORUS_BOX_AY_KEY] = float(ay)
struct.property[Crystal.CHORUS_BOX_AZ_KEY] = float(az)
struct.property[Crystal.CHORUS_BOX_BX_KEY] = float(bx)
struct.property[Crystal.CHORUS_BOX_BY_KEY] = float(by)
struct.property[Crystal.CHORUS_BOX_BZ_KEY] = float(bz)
struct.property[Crystal.CHORUS_BOX_CX_KEY] = float(cx)
struct.property[Crystal.CHORUS_BOX_CY_KEY] = float(cy)
struct.property[Crystal.CHORUS_BOX_CZ_KEY] = float(cz)
[docs]def copy_chorus_box_props(struct1, struct2):
"""
Copy the Chorus box PBC properties from one struct to another if they exist.
If no Chorus properties exist, no error is raised.
:param `structure.Structure` struct1: The structure to copy props from
:param `structure.Structure` struct2: The structure to copy props to
"""
try:
set_pbc_properties(struct2, get_chorus_properties(struct1))
except KeyError:
pass
[docs]def clear_asu_and_fractional_properties(struct):
"""
Clear the atomic ASU and Fractional coordinate properties
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to clear the ASU and fractional properties from
"""
# TODO: This properties are not used anymore anywhere (MATSCI-4673)
F1F2F3_KEY = 's_matsci_f1f2f3'
F1_KEY = 'r_matsci_f1'
F2_KEY = 'r_matsci_f2'
F3_KEY = 'r_matsci_f3'
FRACT_KEYS = [F1_KEY, F2_KEY, F3_KEY]
props = [Crystal.ASU_LABEL_KEY, F1F2F3_KEY] + FRACT_KEYS
for prop in props:
msutils.remove_atom_property(struct, prop)
[docs]def get_cov_radii():
"""
Return a dictionary of atomic covalent radii.
:rtype: dict
:return: dictionary where keys are atomic numbers and
values are atomic covalent radii in Angstrom
"""
global _COV_RADII
if _COV_RADII:
return _COV_RADII
# loop over atoms we will need and build the dictionary
# of covalent radii, do it this way rather than accessing
# by atomic number later so that we can get the max radius
# as well as avoid set up/tear down cycles
mm.mmpdb_initialize(mm.MMERR_DEFAULT_HANDLER)
mmpdb_max_elements = 103
_COV_RADII = {}
for atomic_number in range(1, mmpdb_max_elements + 1):
_COV_RADII[atomic_number] = \
mm.mmpdb_get_cov_radius(atomic_number, mm.MMPDB_VERSION_3)
mm.mmpdb_terminate()
return _COV_RADII
[docs]def get_cov_radius(st, idx, is_cg=None):
"""
Return the covalent radius in Angstrom of the given atom or
coarse grain particle index.
:type st: schrodinger.structure.Structure
:param st: the structure
:type idx: int
:param idx: the atom or coarse grain particle index
:type is_cg: bool or None
:param is_cg: True, if structure is CG, otherwise False. If None, structure
will be evaluated (slow)
:rtype: float
:return: the covalent radius in Angstrom
"""
if is_cg is None:
is_cg = msutils.is_coarse_grain(st)
if is_cg:
# CG particles do not currently have a covalent radius
# property but rather only a VDW radius property, similar
# to coarsegrain.setCGBondLengths just use the VDW radius
return st.atom[idx].radius
else:
atomistic_cov_radii = get_cov_radii()
# allow unknown or dummy atoms by setting the covalent radius
# to 1 (MATSCI-2613)
return atomistic_cov_radii.get(st.atom[idx].atomic_number, 1)
[docs]def get_max_cov_radius(st, is_cg=None):
"""
For atomic input return the maximum covalent radius in Angstrom
over all atoms in the periodic table, for CG input return the same
but over all particles in the given structure.
:type st: schrodinger.structure.Structure
:param st: the structure
:type is_cg: bool or None
:param is_cg: True, if structure is CG, otherwise False. If None, structure
will be evaluated (slow)
:rtype: float
:return: the maximum covalent radius in Angstrom
"""
if is_cg is None:
is_cg = msutils.is_coarse_grain(st)
if is_cg:
return max(get_cov_radius(st, atom.index, is_cg) for atom in st.atom)
else:
atomistic_cov_radii = get_cov_radii()
return max(atomistic_cov_radii.values())
[docs]def get_cov_parameter(st, atomistic_cov_param, is_cg=None):
"""
Return the covalent parameter in Angstrom for the
given structure.
:type st: schrodinger.structure.Structure
:param st: the structure
:type atomistic_cov_param: float
:param atomistic_cov_param: the atomistic covalent parameter
in Angstrom
:type is_cg: bool or None
:param is_cg: True, if structure is CG, otherwise False. If None, structure
will be evaluated (slow)
:rtype: float
:return: the covalent parameter in Angstrom
"""
if is_cg is None:
is_cg = msutils.is_coarse_grain(st)
if not is_cg:
return atomistic_cov_param
else:
# model the corresponding covalent parameter for
# CG systems by using ratios:
#
# atomistic_cov_param / atomistic_ref_cov_radius =
# cg_cov_param / cg_ref_cov_radius
#
# where the reference covalent radii are minimal values, i.e.
# that of hydrogen for atomistic and the default for CG (as
# typically they go up from there), as in get_cov_radius CG
# covalent radii are the same as VDW radii
atomistic_ref_cov_radius = get_cov_radii()[1]
cg_ref_cov_radius = coarsegrain.DEFAULT_RADIUS
return atomistic_cov_param * cg_ref_cov_radius / atomistic_ref_cov_radius
[docs]def get_params_from_vectors(a_vec, b_vec, c_vec):
"""
Return the lattice parameters from the given lattice vectors.
: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: list
:return: contains the a, b, c, alpha, beta, and gamma parameters
"""
pbc = infrastructure.PBC(*a_vec, *b_vec, *c_vec)
return list(pbc.getBoxLengths() + pbc.getBoxAngles())
[docs]def get_lattice_vectors(a_param, b_param, c_param, alpha_param, beta_param,
gamma_param):
"""
Get the lattice vectors of the specified parallelepiped.
:type a_param: float
:param a_param: the length of the parallelepiped along edge a
:type b_param: float
:param b_param: the length of the parallelepiped along edge b
:type c_param: float
:param c_param: the length of the parallelepiped along edge c
:type alpha_param: float
:param alpha_param: the angle between edges b and c
:type beta_param: float
:param beta_param: the angle between edges a and c
:type gamma_param: float
:param gamma_param: the angle between edges a and b
:rtype: numpy.array, numpy.array, numpy.array
:return: the lattice vectors of the parallelepiped
"""
def approx1(value):
return abs(value - 1) < 0.001
# MAE-12963 - Reject NMR structures (a = b = c = 1.0)
if all(map(approx1, (a_param, b_param, c_param))):
raise ValueError(
"All box lengths are near 1. Are these box sizes from NMR?")
# get lattice vectors for unconstrained space group P1
pbc = infrastructure.PBC(a_param, b_param, c_param, alpha_param, beta_param,
gamma_param)
return [numpy.array(v) for v in pbc.getBoxVectors()]
[docs]class ParserWrapper(object):
"""
Manages the argparse module to parse user command line
arguments.
"""
SPACE_GROUP = None
A_PARAM = None
B_PARAM = None
C_PARAM = None
ALPHA_PARAM = None
BETA_PARAM = None
GAMMA_PARAM = None
MIN_ANGLE = 0.0
MAX_ANGLE = 180.0
NCELLA = 1
NCELLB = 1
NCELLC = 1
ORIGIN = [0.0, 0.0, 0.0]
CHOICE_ON = 'on'
CHOICE_OFF = 'off'
BONDING_CHOICES = [CHOICE_ON, CHOICE_OFF]
BONDING_DEFAULT = CHOICE_OFF
BOND_ORDERS_CHOICES = [CHOICE_ON, CHOICE_OFF]
BOND_ORDERS_DEFAULT = CHOICE_OFF
TRANSLATE_CHOICES = list(BONDING_CHOICES)
TRANSLATE_DEFAULT = CHOICE_OFF
COV_MIN = 0.40
COV_OFFSET = 0.40
COV_FACTOR = 1.00
PRINT_INFO = False
PBC_BONDING = True
[docs] def __init__(self, args, scriptname, description):
"""
Create a ParserWrapper object and process it.
:type args: tuple
:param args: command line arguments
:type scriptname: str
:param scriptname: name of this script
:type description: str
:param description: description of this script
"""
# create argparse.ArgumentParser, load it with
# parsables, parse command line args, and return options.
name = '$SCHRODINGER/run ' + scriptname
self.parserobj = argparse.ArgumentParser(
prog=name,
description=description,
add_help=False,
formatter_class=argparse.ArgumentDefaultsHelpFormatter)
self.loadIt()
self.parseArgs(args)
self.validate()
[docs] def loadIt(self):
"""
Load ParserWrapper with options.
"""
space_group_help = """Specify the Hermann-Mauguin symbol
of the space group. The full symbol can be passed or if
the short symbol is unique, for example F m -3 m, that
can be passed as well. If a short symbol is used but
it is non-unique, perhaps because multiple space group
settings are using the same symbol, for example P 2/c, then
the first space group with this short symbol in the list
provided with -print_info will be used. Often this will
be the standard setting of the space group but it is not
guaranteed."""
self.parserobj.add_argument('-space_group',
type=str,
default=self.SPACE_GROUP,
help=space_group_help)
a_help = """Specify the lattice parameter \'a\' in Angstrom.
Lattice lengths should be > 0."""
self.parserobj.add_argument('-a',
type=float,
default=self.A_PARAM,
dest='a_param',
help=a_help)
b_help = """Specify the lattice parameter \'b\' in Angstrom.
Lattice lengths should be > 0."""
self.parserobj.add_argument('-b',
type=float,
default=self.B_PARAM,
dest='b_param',
help=b_help)
c_help = """Specify the lattice parameter \'c\' in Angstrom.
Lattice lengths should be > 0."""
self.parserobj.add_argument('-c',
type=float,
default=self.C_PARAM,
dest='c_param',
help=c_help)
alpha_help = """Specify the lattice parameter \'alpha\' in
degrees. Lattice angles should be in (%s, %s)."""
self.parserobj.add_argument('-alpha', type=float,
default=self.ALPHA_PARAM, dest='alpha_param', help=alpha_help % \
(self.MIN_ANGLE, self.MAX_ANGLE))
beta_help = """Specify the lattice parameter \'beta\' in
degrees. Lattice angles should be in (%s, %s)."""
self.parserobj.add_argument('-beta', type=float,
default=self.BETA_PARAM, dest='beta_param', help=beta_help % \
(self.MIN_ANGLE, self.MAX_ANGLE))
gamma_help = """Specify the lattice parameter \'gamma\' in
degrees. Lattice angles should be in (%s, %s)."""
self.parserobj.add_argument('-gamma', type=float,
default=self.GAMMA_PARAM, dest='gamma_param', help=gamma_help % \
(self.MIN_ANGLE, self.MAX_ANGLE))
ncella_help = """Specify the desired number of unit cells along
lattice vector \'a\'."""
self.parserobj.add_argument('-ncella',
type=int,
default=self.NCELLA,
help=ncella_help)
ncellb_help = """Specify the desired number of unit cells along
lattice vector \'b\'."""
self.parserobj.add_argument('-ncellb',
type=int,
default=self.NCELLB,
help=ncellb_help)
ncellc_help = """Specify the desired number of unit cells along
lattice vector \'c\'."""
self.parserobj.add_argument('-ncellc',
type=int,
default=self.NCELLC,
help=ncellc_help)
origin_help = """By default the origin of the unit cell will be
(0, 0, 0). If a translation to the first unit cell is being
performed then use this option to define where that unit cell
is. Note that this option is not merely translating a given
unit cell but rather defining the reference frame. It
effectively allows one to cut different valid unit cells from the
larger periodic structure. This option is not used if a
translation is not being performed. Specify it in fractional
coordinates. The 'anchor_x_y_z' (x, y, and z in Cartesian
coordinates) value of the 's_mae_pbc_position' structure
property on the input ASU can also be used to specify an origin
though this origin flag will take precedence."""
self.parserobj.add_argument('-origin',
type=float,
nargs=3,
metavar='FLOAT',
help=origin_help)
bonding_help = """Specify if this script should handle the bonding both
within the unit cell and between unit cells. See the option -bond_orders
for information about assigning bond orders."""
self.parserobj.add_argument('-bonding',
type=str,
choices=self.BONDING_CHOICES,
default=self.BONDING_DEFAULT,
help=bonding_help)
bond_orders_help = """Specify if this script should assign bond orders
both within the unit cell and between unit cells."""
self.parserobj.add_argument('-bond_orders',
type=str,
choices=self.BOND_ORDERS_CHOICES,
default=self.BOND_ORDERS_DEFAULT,
help=bond_orders_help)
translate_help = """Specify if this script should translate all atoms
to the first unit cell which can be defined using the -origin
option."""
self.parserobj.add_argument('-translate',
type=str,
choices=self.TRANSLATE_CHOICES,
default=self.TRANSLATE_DEFAULT,
help=translate_help)
translatec_help = """Specify if this script should translate all
molecules to the first unit cell. Cannot be used together with
-translate"""
self.parserobj.add_argument('-translate_centroids',
action='store_true',
help=translatec_help)
cov_offset_help = """Bonds will be drawn using the covalent radii
of atoms. The maximum distance for a connection is the sum of
the radii of the two atoms plus this offset in Angstrom.
Increasing this value will increase the number of bonds. Note
that the covalent radii used will be returned as atom properties
of the structure."""
self.parserobj.add_argument('-cov_offset',
type=float,
default=self.COV_OFFSET,
help=cov_offset_help)
fract_offset_help = """The following threshold is used when making
comparisions involving floating point fractional coordinate
values. Many tests indicate that this is an acceptable value
however if your system seems to have less-than or more-than
the correct number of atoms, for example if the stoichiometry
is incorrect, then it is possible that adjusting this value
will correct this behavior."""
self.parserobj.add_argument('-fract_offset',
type=float,
default=FRACT_OFFSET,
help=fract_offset_help)
print_info_help = """Print detailed information about all space
groups to a log file and exit. For example use this option
to see a list of supported short and full Hermann-Mauguin
symbols. If you want to also log the symmetry operators of
the space groups then in addition pass the -verbose flag. They
will then be at the bottom of the log file."""
self.parserobj.add_argument('-print_info',
action='store_true',
help=print_info_help)
LATTICE_PARAM_KEYS = [
Crystal.SPACE_GROUP_KEY, Crystal.A_KEY, Crystal.B_KEY,
Crystal.C_KEY, Crystal.ALPHA_KEY, Crystal.BETA_KEY,
Crystal.GAMMA_KEY
]
input_file_help = """Specify a file containing the asymmetric unit of
the crystal (Maestro, CIF, PDB). Some script parameters may be
specified in this input file, rather than by the script arguments
given below, by providing the following list of structure
properties: %s. Note that script arguments will always take
precedence.""" % ', '.join(LATTICE_PARAM_KEYS)
self.parserobj.add_argument('input_file',
type=str,
nargs='?',
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
"""
# move input files, which will be the only positional arguments,
# to be first in the list of args
args = [str(arg) for arg in args]
infiles = []
for arg in list(args):
if fileutils.is_maestro_file(arg):
args.remove(arg)
infiles.append(arg)
args = infiles + args
self.options = self.parserobj.parse_args(args)
[docs] def validate(self):
"""Validate parser options."""
if (self.options.translate == self.CHOICE_ON and
self.options.translate_centroids):
self.parserobj.error('Flags -translate and -translate_centroids '
'cannot be used together.')
[docs]def get_duplicate_atoms(struct,
pbc=None,
atoms_to_check=None,
duplicate_thresh=OVERLAP_ATOM_THRESHOLD):
"""
Get atoms to keep and atoms to delete from a structure that possibly
contains duplicate atoms that are within the defined threshold. Of the
redundant atoms that with the lowest index is kept.
:type struct: schrodinger.Structure.structure
:param struct: the structure object to find duplicate atoms
:type pbc: `schrodinger.infra.structure.PBC` or None
:param pbc: The infrastructure PBC created from the Chorus box properties,
if None, system will be treated as finite
:type atoms_to_check: set
:param atoms_to_check: indices of atoms that are to be checked for
duplicate copies
:type duplicate_thresh: float
:param duplicate_thresh: distance used to define duplicate atoms,
conservatively chosen based on considering the resolutions of some
crystal structures and the wavelengths of light used to obtain the
diffraction patterns
:rtype: list, list
:return: First list of atom indices to keep. Second is list of lists of atom
indices to be removed. Each item in the list is a list of duplicated
atom indices corresponding to the atom index from to_keep list
"""
if atoms_to_check is None:
atoms_to_check = set(range(1, struct.atom_total + 1))
if pbc:
cell = infrastructure.DistanceCell(struct, duplicate_thresh, pbc)
else:
cell = infrastructure.DistanceCell(struct, duplicate_thresh)
# next from the list of duplicates for each requested atom
# keep the one with the lowest atomic index and mark the rest
# for deletion
atoms_to_delete = []
atoms_to_keep = []
while atoms_to_check:
idx = atoms_to_check.pop()
atom = struct.atom[idx]
matches = cell.query_atoms(*atom.xyz)
indices = sorted(match.getIndex() for match in matches)
if len(indices) < 2:
continue
atoms_to_keep.append(indices[0])
atoms_to_delete.append(indices[1:])
atoms_to_check.difference_update(set(indices[1:]))
return atoms_to_keep, atoms_to_delete
[docs]def delete_duplicate_atoms(astructure,
atoms_to_check=None,
duplicate_thresh=OVERLAP_ATOM_THRESHOLD,
fract_offset=FRACT_OFFSET,
preserve_bonding=False):
"""
Delete duplicate atoms that are within the defined threshold. Of
the redundant atoms that with the lowest index is kept. If transform
is not None then this function will use the periodic boundary conditions
defined in transform when determining redundant atoms.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object from which the duplicate
atoms will be deleted
:type atoms_to_check: set
:param atoms_to_check: indices of atoms that are to be checked for
duplicate copies, for a given redundant set of atoms any index can
be passed in
:type duplicate_thresh: float
:param duplicate_thresh: distance used to define duplicate atoms,
conservatively chosen based on considering the resolutions of some
crystal structures and the wavelengths of light used to obtain the
diffraction patterns
:type fract_offset: float
:param fract_offset: the threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:type preserve_bonding: bool
:param preserve_bonding: If True, preserve bonding between atoms
(might be slow)
:rtype: dict
:return: renumber_map, keys are original indices, values are new
indices or None if the atom was deleted
"""
tmpst = astructure.copy()
# finalize the list of atoms for which duplicates will be checked
if not atoms_to_check:
atoms_to_check = set(range(1, tmpst.atom_total + 1))
# move atoms into the cell
translate_to_cell(tmpst, fract_offset=fract_offset)
try:
# Prioritize lattice parameters over chorus
params = get_lattice_param_properties(tmpst)
pbc = infrastructure.PBC(*get_chorus_from_params(params))
except KeyError:
try:
pbc = infrastructure.PBC(*get_chorus_properties(tmpst))
except KeyError:
pbc = None
if pbc:
cell = infrastructure.DistanceCell(tmpst, duplicate_thresh, pbc)
else:
cell = infrastructure.DistanceCell(tmpst, duplicate_thresh)
atoms_to_keep, atoms_to_delete = get_duplicate_atoms(
tmpst,
pbc=pbc,
atoms_to_check=atoms_to_check,
duplicate_thresh=duplicate_thresh)
renumber_map = dict()
if len(atoms_to_delete):
if preserve_bonding:
preserve_bonds(astructure, atoms_to_keep, atoms_to_delete)
renumber_map = astructure.deleteAtoms(set(
numpy.concatenate(atoms_to_delete)),
renumber_map=True)
return renumber_map
[docs]def max_connect_distance(cov_rad_a,
cov_rad_b,
cov_factor=ParserWrapper.COV_FACTOR,
cov_offset=ParserWrapper.COV_OFFSET):
"""
Return the maximum bonding distance for the given covalent
radii and distance equation parameters.
:type cov_rad_a: float
:param cov_rad_a: covalent radii for atom a
:type cov_rad_b: float
:param cov_rad_b: covalent radii for atom b
:type cov_factor: float
:param cov_factor: the maximum distance for a connection is
the sum of the covalent radii of the two atoms weighted by this factor
plus the cov_offset in angstrom, increasing this value will increase
the number of connections, this value is unit-less
:type cov_offset: float
:param cov_offset: the maximum distance for a connection is
the sum of the covalent radii of the two atoms weighted by cov_factor
plus this offset in angstrom, increasing this value will increase
the number of connections
:rtype: float
:return: the maximum bonding distance
"""
# depending on cov_offset and cov_factor the bonding for
# materials problems can be poor. increasing or decreasing
# these values can provide the desired result but only at the
# expense of sacrificing some other result which may not be
# desired. i think we will have to control this behavior with
# either a modified bonding protocol for materials or by
# accumulating a database of radii for different types of
# atoms with different oxidation states, etc., the maxium
# bonding distance, y, used herein is general, i.e.
#
# y = cov_factor*(cov_radii_a + cov_radii_b) + cov_offset
#
# the reason for this is because our existing connectors use
# either one or the other, for example OpenBabel will use the
# cov_offset while mmjag and mmpdb both use cov_factor, for
# materials it looks like cov_offset is the way to do it but
# that deserves more testing
return cov_factor * (cov_rad_a + cov_rad_b) + cov_offset
[docs]def connect_atoms(astructure,
atoms_to_connect=None,
cov_min=ParserWrapper.COV_MIN,
cov_offset=ParserWrapper.COV_OFFSET,
cov_factor=ParserWrapper.COV_FACTOR,
delete_existing=True,
cov_radii_props=False,
pbc_bonding=ParserWrapper.PBC_BONDING,
only_pbc_bonds=False,
check_also_reg_bond=None,
max_valencies=None):
"""
Connect the atoms in a structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object for which the connections
are wanted
:type atoms_to_connect: list
:param atoms_to_connect: atom indices to consider connecting, if
None then all atoms will be considered
:type cov_min: float
:param cov_min: the minimum distance for a connection in angstrom
:type cov_offset: float
:param cov_offset: the maximum distance for a connection is
the sum of the covalent radii of the two atoms weighted by cov_factor
plus this offset in angstrom, increasing this value will increase
the number of connections
:type cov_factor: float
:param cov_factor: the maximum distance for a connection is
the sum of the covalent radii of the two atoms weighted by this factor
plus the cov_offset in angstrom, increasing this value will increase
the number of connections, this value is unit-less
:type delete_existing: bool
:param delete_existing: indicates whether existing connections
should first be deleted, i.e. any bonds to atoms in
atoms_to_connect will be deleted
:type cov_radii_props: bool
:param cov_radii_props: set the atomic covalent radii used in the
connectivity protocol as atom properties. Enabling this, causes a
decrease in performance of approx. 10%.
:type pbc_bonding: bool
:param pbc_bonding: if True and if the input structure has the
chorus properties defined, then PBC bonds will be created
:type only_pbc_bonds: bool
:param only_pbc_bonds: if True then only PBC bonds will be
created, i.e. regular bonds will not be created
:type check_also_reg_bond: bool or None
:param check_also_reg_bond: if True then PBC bonds will be
checked to see if they are also regular bonds, meaning that
one of the atoms of the PBC bond is covalently bound to two
copies of the other atom, one inside the cell and one outside
the cell. If None, this will be evaluated
:type max_valencies: dict or None
:param max_valencies: The maximum valencies to use for each element.
If None, the valencies defined in get_max_valencies() will be used.
:rtype: list, dict, dict
:return: maximally_bonded_atoms, contains the indices of any
maximally bonded atoms, and two dictionaries (one for regular
bonds and one for PBC bonds) where keys are tuples of bonding
atom indices and values are bond orders (for regular bonds) or
tuples of bond orders and booleans indicating whether the PBC
bond is also a regular bond (for PBC bonds)
"""
if max_valencies is None:
max_valencies = elementalprops.get_max_valencies()
if only_pbc_bonds and not pbc_bonding:
return [], {}, {}
is_cg = msutils.is_coarse_grain(astructure)
cov_min = get_cov_parameter(astructure, cov_min, is_cg=is_cg)
cov_offset = get_cov_parameter(astructure, cov_offset, is_cg=is_cg)
# get the covalent radius and optionally set it as a property
def handle_radii(atom):
cov_rad = get_cov_radius(astructure, atom.index, is_cg=is_cg)
if cov_radii_props:
atom.property[COV_RADIUS_KEY] = cov_rad
return cov_rad
def is_maximally_bonded(atom):
if atom.bond_total == max_valencies[atom.element]:
maximally_bonded_atoms.add(atom.index)
return True
# prepare, see MATSCI-5996 where it was decided to honor incoming
# bond properties
if not atoms_to_connect:
atoms_to_connect = list(range(1, astructure.atom_total + 1))
bonds_w_props = {}
if delete_existing:
for bond in astructure.bond:
idxs = (bond.atom1.index, bond.atom2.index)
bonds_w_props[idxs] = list(bond.property.items())
build.delete_bonds_to_atoms(
(astructure.atom[idx] for idx in atoms_to_connect))
# get chorus properties
try:
chorus_properties = get_chorus_properties(astructure)
pbc = infrastructure.PBC(*chorus_properties)
except KeyError:
chorus_properties = None
pbc = None
# determine check_also_reg_bond
if check_also_reg_bond is None:
if chorus_properties and pbc_bonding:
check_also_reg_bond = get_check_also_reg_bond(astructure,
pbc=pbc,
is_cg=is_cg)
else:
check_also_reg_bond = False
# the lattice vectors will be needed to check if PBC bonds
# are also regular bonds
if chorus_properties and pbc_bonding and check_also_reg_bond:
unit_lattice_vectors = get_unit_lattice_vector_info(astructure)
else:
unit_lattice_vectors = None
# get potential bonding pairs for the requested atoms using
# a distance cell, for the distance use the max distance
# equation given the atom with the largest covalent radius
if len(atoms_to_connect) > 2:
max_cov_radius = get_max_cov_radius(astructure, is_cg=is_cg)
cell_distance = max_connect_distance(max_cov_radius,
max_cov_radius,
cov_factor=cov_factor,
cov_offset=cov_offset)
pairs = get_cell_pairs(astructure,
cell_distance,
pbc_bonding=pbc_bonding,
atom_indices=atoms_to_connect,
chorus_properties=chorus_properties)
else:
pairs = set()
if len(atoms_to_connect) == 2:
pairs.add(tuple(atoms_to_connect))
# do the bonding over the pairs within cell_distance
maximally_bonded_atoms = set()
bonds_made = {}
pbc_bonds_made = {}
for pair in pairs:
a_atom, b_atom = [astructure.atom[x] for x in pair]
if (a_atom in maximally_bonded_atoms or
b_atom in maximally_bonded_atoms or
is_maximally_bonded(a_atom) or is_maximally_bonded(b_atom)):
continue
# get the distance
pbc_bond = False
also_reg_bond = False
if chorus_properties and pbc_bonding:
pbc_bond, also_reg_bond, distance = \
is_pbc_bond(astructure, a_atom, b_atom, \
check_also_reg_bond=check_also_reg_bond, \
unit_lattice_vectors=unit_lattice_vectors,
pbc=pbc)
else:
distance = mm.mmct_atom_get_distance(astructure.handle, a_atom,
astructure.handle, b_atom)
if not pbc_bond and only_pbc_bonds:
continue
# compare against the min and max distances
cov_rad_a, cov_rad_b = list(map(handle_radii, [a_atom, b_atom]))
max_distance = max_connect_distance(cov_rad_a, cov_rad_b, \
cov_factor=cov_factor, cov_offset=cov_offset)
if cov_min < distance < max_distance:
bond_order = 1
add_labeled_pbc_bond(astructure,
a_atom,
b_atom,
bond_order,
is_pbc_bond=pbc_bond,
also_reg_bond=also_reg_bond,
color=None)
if pbc_bond:
pbc_bonds_made[pair] = (bond_order, also_reg_bond)
else:
bonds_made[pair] = bond_order
reserved_pbc_bond_keys = (PBCBond.PBC_BOND_KEY, PBCBond.ALSO_REG_BOND_KEY,
PBCBond.PBC_BOND_COLOR_KEY)
for idxs, props in bonds_w_props.items():
bond = astructure.getBond(*idxs)
if bond:
for key, value in props:
# in the interest of honoring only true meta-data
# as opposed to infrastructure that should not be
# changed skip Maestro properties, most importantly
# 'i_m_order' which will depend if assigning bond
# orders in this script among other things, also
# bond representation 'i_m_from_rep' and 'i_m_to_rep'
# which are set above, lastly skip PBC bond properties
# which are also set above
if key[1:].startswith('_m_') or key in reserved_pbc_bond_keys:
continue
bond.property[key] = value
maximally_bonded_atoms = list(maximally_bonded_atoms)
maximally_bonded_atoms.sort()
return maximally_bonded_atoms, bonds_made, pbc_bonds_made
# the following five functions are used to classify fractional
# coordinate values as (1) before lower edge, (2) on lower edge,
# (3) inside cell, (4) on upper edge, and (5) after upper edge
[docs]def before_lower_edge(coord, lower_bound):
"""
Return True if this fractional coordinate value is before the
lower edge of the cell.
:type coord: float
:param coord: the fractional coordinate being tested
:type lower_bound: float
:param lower_bound: the lower bound for this fractional
coordinate
:rtype: bool
:return: True if the coordinate is before the lower edge,
False otherwise
"""
if coord < lower_bound:
return True
return False
[docs]def on_lower_edge(coord, lower_bound, fract_offset=FRACT_OFFSET):
"""
Return True if this fractional coordinate value is on the
lower edge of the cell.
:type coord: float
:param coord: the fractional coordinate being tested
:type lower_bound: float
:param lower_bound: the lower bound for this fractional
coordinate
:type fract_offset: float
:param fract_offset: used to make floating point comparisons
:rtype: bool
:return: True if the coordinate is on the lower edge,
False otherwise
"""
if abs(coord - lower_bound) <= fract_offset:
return True
return False
[docs]def inside_cell(coord, lower_bound, upper_bound, fract_offset=FRACT_OFFSET):
"""
Return True if this fractional coordinate value is inside the
cell.
:type coord: float
:param coord: the fractional coordinate being tested
:type lower_bound: float
:param lower_bound: the lower bound for this fractional
coordinate
:type upper_bound: float
:param upper_bound: the upper bound for this fractional
coordinate
:type fract_offset: float
:param fract_offset: used to make floating point comparisons
:rtype: bool
:return: True if the coordinate is inside the cell,
False otherwise
"""
if lower_bound <= coord < upper_bound - fract_offset:
return True
return False
[docs]def on_upper_edge(coord, lower_bound, fract_offset=FRACT_OFFSET):
"""
Return True if this fractional coordinate value is on the
upper edge of the cell.
:type coord: float
:param coord: the fractional coordinate being tested
:type lower_bound: float
:param lower_bound: the lower bound for this fractional
coordinate
:type fract_offset: float
:param fract_offset: used to make floating point comparisons
:rtype: bool
:return: True if the coordinate is on the upper edge,
False otherwise
"""
if abs(coord - lower_bound - 1) <= fract_offset:
return True
return False
[docs]def after_upper_edge(coord, lower_bound, fract_offset=FRACT_OFFSET):
"""
Return True if this fractional coordinate value is after the
upper edge of the cell.
:type coord: float
:param coord: the fractional coordinate being tested
:type lower_bound: float
:param lower_bound: the lower bound for this fractional
coordinate
:type fract_offset: float
:param fract_offset: used to make floating point comparisons
:rtype: bool
:return: True if the coordinate is after the upper edge,
False otherwise
"""
if lower_bound + 1 - fract_offset <= coord:
return True
return False
[docs]def translate_to_cell(astructure,
fract_offset=FRACT_OFFSET,
origin=ParserWrapper.ORIGIN,
extents=None):
"""
Translate the fractional coordinate definitions of the
atoms of the given structure so that they are all in the
cell defined with the given origin and given extents and
optionally also actually transform the Cartesian coordinates
of the atoms using the specified fractional-to-Cartesian
transform.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object whose atoms will be
translated
:type fract_offset: float
:param fract_offset: the threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:type origin: list
:param origin: the origin of the cell to which to translate in
fractional coordinates
:type extents: list
:param extents: the number of cells along each lattice vector
"""
if extents is None:
extents = [1, 1, 1]
f1_lower_bound, f2_lower_bound, f3_lower_bound = origin
f1_upper_bound = f1_lower_bound + extents[0]
f2_upper_bound = f2_lower_bound + extents[1]
f3_upper_bound = f3_lower_bound + extents[2]
# the unit cell is defined on 0 <= n < N, with N = 1, for
# extented cells N > 1
def process_atom(coord, lower_bound, upper_bound):
while before_lower_edge(coord, lower_bound):
coord += upper_bound - lower_bound
while after_upper_edge(coord, upper_bound - 1, fract_offset):
coord -= upper_bound - lower_bound
return coord
vecs = numpy.array(
get_lattice_vectors(*get_lattice_param_properties(astructure)))
fracs = trans_cart_to_frac_from_vecs(astructure.getXYZ(), *vecs)
new_fracs = numpy.array(fracs)
for idx, atom_fracs in enumerate(fracs):
f1 = process_atom(atom_fracs[0], f1_lower_bound, f1_upper_bound)
f2 = process_atom(atom_fracs[1], f2_lower_bound, f2_upper_bound)
f3 = process_atom(atom_fracs[2], f3_lower_bound, f3_upper_bound)
new_fracs[idx] = [f1, f2, f3]
astructure.setXYZ(trans_frac_to_cart_from_vecs(new_fracs, *vecs))
[docs]def get_gcd_list_ints(list_of_ints):
"""
Return the greatest common divisor (GCD) of a list of integers.
:type list_of_ints: list
:param list_of_ints: the list of ints for which the GCD is desired
:rtype: int
:return: gcd, the GCD
"""
agcd = list_of_ints[0]
for an_int in list_of_ints[1:]:
agcd = gcd(agcd, an_int)
return agcd
[docs]def get_carts_from_anchor_string(anchor):
"""
Return the Cartesian coordinates from the given
's_mae_pbc_position' anchor string, i.e. 'anchor_x_y_z'.
:type anchor: str
:param anchor: the anchor string
:raise ValueError: if the given string is not an anchor string
:rtype: list
:return: the Cartesians floats
"""
if not anchor.startswith(ANCHOR_PREFIX):
msg = ('The given string does not start with %s.' % ANCHOR_PREFIX)
raise ValueError(msg)
return [float(i) for i in anchor.split('_')[1:]]
[docs]class Crystal(object):
"""
Main class for generating crystals.
"""
SPACE_GROUP_KEY = SPACE_GROUP_KEY
A_KEY = 'r_pdb_PDB_CRYST1_a'
B_KEY = 'r_pdb_PDB_CRYST1_b'
C_KEY = 'r_pdb_PDB_CRYST1_c'
ALPHA_KEY = 'r_pdb_PDB_CRYST1_alpha'
BETA_KEY = 'r_pdb_PDB_CRYST1_beta'
GAMMA_KEY = 'r_pdb_PDB_CRYST1_gamma'
SPACE_GROUP_ID_KEY = SPACE_GROUP_ID_KEY
CHORUS_BOX_BASE_KEY = 'r_chorus_box_'
CHORUS_BOX_AX_KEY = CHORUS_BOX_BASE_KEY + 'ax'
CHORUS_BOX_AY_KEY = CHORUS_BOX_BASE_KEY + 'ay'
CHORUS_BOX_AZ_KEY = CHORUS_BOX_BASE_KEY + 'az'
CHORUS_BOX_BX_KEY = CHORUS_BOX_BASE_KEY + 'bx'
CHORUS_BOX_BY_KEY = CHORUS_BOX_BASE_KEY + 'by'
CHORUS_BOX_BZ_KEY = CHORUS_BOX_BASE_KEY + 'bz'
CHORUS_BOX_CX_KEY = CHORUS_BOX_BASE_KEY + 'cx'
CHORUS_BOX_CY_KEY = CHORUS_BOX_BASE_KEY + 'cy'
CHORUS_BOX_CZ_KEY = CHORUS_BOX_BASE_KEY + 'cz'
CHORUS_BOX_A_KEYS = [
CHORUS_BOX_AX_KEY, CHORUS_BOX_AY_KEY, CHORUS_BOX_AZ_KEY
]
CHORUS_BOX_B_KEYS = [
CHORUS_BOX_BX_KEY, CHORUS_BOX_BY_KEY, CHORUS_BOX_BZ_KEY
]
CHORUS_BOX_C_KEYS = [
CHORUS_BOX_CX_KEY, CHORUS_BOX_CY_KEY, CHORUS_BOX_CZ_KEY
]
# keys is all one list, key_vectors is a list of lists
CHORUS_BOX_KEYS = CHORUS_BOX_A_KEYS + CHORUS_BOX_B_KEYS + CHORUS_BOX_C_KEYS
CHORUS_BOX_KEY_VECTORS = [
CHORUS_BOX_A_KEYS, CHORUS_BOX_B_KEYS, CHORUS_BOX_C_KEYS
]
NUM_DECIMAL_COORDS = 4
MSGWIDTH = 80
# basis column vectors
GENERAL_VEC_1 = numpy.array([1.0, 0.0, 0.0])
GENERAL_VEC_2 = numpy.array([0.0, 1.0, 0.0])
GENERAL_VEC_3 = numpy.array([0.0, 0.0, 1.0])
SYMMETRY_LABEL_KEY = 'i_matsci_Symmetry_Label'
ASU_LABEL_KEY = 'b_matsci_ASU_Atom'
CM_TO_ANGSTROM = float(100000000)
# making these Maestro properties (see MATSCI-1216)
UNIT_CELL_FORMULA_KEY = msprops.UNIT_CELL_FORMULA_KEY
UNIT_CELL_VOLUME_KEY = msprops.UNIT_CELL_VOLUME_KEY
UNIT_CELL_DENSITY_KEY = msprops.UNIT_CELL_DENSITY_KEY
[docs] def __init__(self,
asymmetric_unit,
space_group=ParserWrapper.SPACE_GROUP,
a_param=ParserWrapper.A_PARAM,
b_param=ParserWrapper.B_PARAM,
c_param=ParserWrapper.C_PARAM,
alpha_param=ParserWrapper.ALPHA_PARAM,
beta_param=ParserWrapper.BETA_PARAM,
gamma_param=ParserWrapper.GAMMA_PARAM,
ncella=ParserWrapper.NCELLA,
ncellb=ParserWrapper.NCELLB,
ncellc=ParserWrapper.NCELLC,
origin=None,
bonding=ParserWrapper.BONDING_DEFAULT,
bond_orders=ParserWrapper.BOND_ORDERS_DEFAULT,
translate=ParserWrapper.TRANSLATE_DEFAULT,
translate_centroids=False,
cov_offset=ParserWrapper.COV_OFFSET,
fract_offset=FRACT_OFFSET,
overlap_tresh=OVERLAP_ATOM_THRESHOLD,
use_existing_pbc_bonds=False,
logger=None,
valency_exceptions=None):
"""
Create an instance.
:type asymmetric_unit: schrodinger.Structure.structure
:param asymmetric_unit: the ASU
:type space_group: str
:param space_group: the full or short Hermann-Mauguin symbol
of the space group
:type a_param: float
:param a_param: the lattice parameter a in Angstrom
:type b_param: float
:param b_param: the lattice parameter b in Angstrom
:type c_param: float
:param c_param: the lattice parameter c in Angstrom
:type alpha_param: float
:param alpha_param: the lattice parameter alpha in degrees
:type beta_param: float
:param beta_param: the lattice parameter beta in degrees
:type gamma_param: float
:param gamma_param: the lattice parameter gamma in degrees
:type ncella: int
:param ncella: the number of unit cells to generate
along lattice vector a
:type ncellb: int
:param ncellb: the number of unit cells to generate
along lattice vector b
:type ncellc: int
:param ncellc: the number of unit cells to generate
along lattice vector c
:type origin: list or None
:param origin: the origin of the unit cell, i.e. list of three floats,
in fractional coordinates or None if one is not provided
:type bonding: str
:param bonding: specifies how the bonding should be handled, takes
a string in BONDING_CHOICES
:type bond_orders: str
:param bond_orders: specifies how the bond orders should be handled,
takes a string in BOND_ORDERS_CHOICES
:type translate: str
:param translate: specifies how the translation to the first unit cell
should be handled, takes a string in TRANSLATE_CHOICES
:param bool translate_centroids: If True, after the cell is built,
entire molecules are moved into the first cell. Cannot be used
together with translate
:type cov_offset: float
:param cov_offset: the maximum distance for drawn bonds is the sum
of the covalent radii of the two atoms plus this offset in angstrom
:type fract_offset: float
:param fract_offset: the threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:param float overlap_tresh: distance used to define overlapping atoms.
If 0, no atoms will be deleted
:type use_existing_pbc_bonds: bool
:param use_existing_pbc_bonds: rather than recalculating PBC bonds
use existing PBC bonds
:type logger: logging.getLogger
:param logger: output logger
:param dict valency_exceptions: The exceptions to update default
valencies with.
"""
# Both cannot be True
assert not (translate == ParserWrapper.CHOICE_ON and
translate_centroids)
self.asymmetric_unit = asymmetric_unit
self.space_group = space_group
self.a_param = a_param
self.b_param = b_param
self.c_param = c_param
self.alpha_param = alpha_param
self.beta_param = beta_param
self.gamma_param = gamma_param
self.ncella = ncella
self.ncellb = ncellb
self.ncellc = ncellc
self.origin = origin
self.bonding = bonding
self.bond_orders = bond_orders
self.translate = translate
self.translate_centroids = translate_centroids
self.cov_offset = cov_offset
self.fract_offset = fract_offset
self.overlap_tresh = overlap_tresh
self.use_existing_pbc_bonds = use_existing_pbc_bonds
self.logger = logger
self.do_bonding = False
self.do_bond_orders = False
self.do_translation = False
self.xyz2cry = None
self.cry2xyz = None
self.a_latt_vec_cart_basis = None
self.b_latt_vec_cart_basis = None
self.c_latt_vec_cart_basis = None
self.ra_latt_vec_cart_basis = None
self.rb_latt_vec_cart_basis = None
self.rc_latt_vec_cart_basis = None
self.spgobj = None
self.crystal_unit_cell = None
self.num_sym_unique = None
self.unit_cell_formula = None
self.unit_cell_volume = None
self.unit_cell_density = None
self.pbc_bonds_by_cell = {}
self.indicies_by_cell = {}
self.crystal_super_cell = None
# Copy, since it is updated
self.max_valencies = dict(elementalprops.get_max_valencies())
if valency_exceptions:
self.max_valencies.update(valency_exceptions)
[docs] def handleSettings(self):
"""
Handle the settings of this crystal build.
"""
# determine if this script should handle the bonding
self.do_bonding = self.doBonding()
# determine if a translation to the first unit cell should be performed
self.do_translation = self.doTranslation()
# determine if we are assigning bond orders
self.do_bond_orders = self.doBondOrders()
[docs] def doBonding(self):
"""
Determine if this script should handle the bonding
both inside and between unit cells.
:rtype: bool
:return: True if this script should handle the bonding,
False otherwise.
"""
return self.bonding == ParserWrapper.CHOICE_ON
[docs] def doBondOrders(self):
"""
Determine if this script should handle the assignment
of bond orders both inside and between unit cells.
:rtype: bool
:return: True if this script should handle the assignment,
False otherwise.
"""
return self.bond_orders == ParserWrapper.CHOICE_ON
[docs] def doTranslation(self):
"""
Determine if a translation to the first unit cell should be performed.
:rtype: bool
:return: True if this translation should be performed,
False otherwise.
"""
return self.translate == ParserWrapper.CHOICE_ON
[docs] def getLatticeParameters(self):
"""
Get lattice parameters.
"""
def get_prop(key):
return self.asymmetric_unit.property.get(key)
# some properties can be obtained from structure level
# properties rather than from command line arguments
if self.space_group is ParserWrapper.SPACE_GROUP:
self.space_group = get_prop(self.SPACE_GROUP_KEY)
if self.a_param is ParserWrapper.A_PARAM:
self.a_param = get_prop(self.A_KEY)
if self.b_param is ParserWrapper.B_PARAM:
self.b_param = get_prop(self.B_KEY)
if self.c_param is ParserWrapper.C_PARAM:
self.c_param = get_prop(self.C_KEY)
if self.alpha_param is ParserWrapper.ALPHA_PARAM:
self.alpha_param = get_prop(self.ALPHA_KEY)
if self.beta_param is ParserWrapper.BETA_PARAM:
self.beta_param = get_prop(self.BETA_KEY)
if self.gamma_param is ParserWrapper.GAMMA_PARAM:
self.gamma_param = get_prop(self.GAMMA_KEY)
[docs] def updateLatticeProperties(self):
"""
Update the lattice properties for the asymmetric unit
structure.
"""
self.asymmetric_unit.property[self.SPACE_GROUP_KEY] = self.space_group
lattice_properties = [self.a_param, self.b_param, self.c_param, \
self.alpha_param, self.beta_param, self.gamma_param]
set_lattice_properties(self.asymmetric_unit, lattice_properties)
def _setCrystalSymmetryFromHandle(self, mmcrystal_handle):
"""
Set the crystal symmetry from the given handle.
:type mmcrystal_handle: handle
:param mmcrystal_handle: the mmcrystal handle
"""
# the important stuff is in the cryxyz and xyz2cry (fractional
# to Cartesian coordinates and back) transforms
nsym, self.cry2xyz, self.xyz2cry = \
mm.mmcrystal_set_symmetry(mmcrystal_handle, self.space_group,
self.a_param, self.b_param, self.c_param,
self.alpha_param, self.beta_param, self.gamma_param)
[docs] def setCrystalSymmetry(self):
"""
Set the crystal symmetry.
"""
mmcrystal_handle = mm.mmcrystal_new()
try:
self._setCrystalSymmetryFromHandle(mmcrystal_handle)
finally:
mm.mmcrystal_delete(mmcrystal_handle)
[docs] def determineBasisVectors(self):
"""
Determine the lattice vectors in the cartesian basis and the
cartesian vectors in the lattice basis. Also determine the
reciprocal lattice vectors.
"""
# update transforms
self.cry2xyz = numpy.array(self.cry2xyz)
self.xyz2cry = numpy.array(self.xyz2cry)
# get bases
a_latt_vec_latt_basis = i_cart_vec_cart_basis = self.GENERAL_VEC_1
b_latt_vec_latt_basis = j_cart_vec_cart_basis = self.GENERAL_VEC_2
c_latt_vec_latt_basis = k_cart_vec_cart_basis = self.GENERAL_VEC_3
# get lattice vectors in cartesian basis, these are not symmetric
# matricies and so we must perform the transform from the left on
# a column vector or apply the transpose from the right on a row
# vector
self.a_latt_vec_cart_basis = numpy.dot(self.cry2xyz,
a_latt_vec_latt_basis)
self.b_latt_vec_cart_basis = numpy.dot(self.cry2xyz,
b_latt_vec_latt_basis)
self.c_latt_vec_cart_basis = numpy.dot(self.cry2xyz,
c_latt_vec_latt_basis)
# get reciprocal lattice vectors, use the same column vector conventions
ra_vec, rb_vec, rc_vec = get_reciprocal_lattice_vectors(self.a_latt_vec_cart_basis, \
self.b_latt_vec_cart_basis, self.c_latt_vec_cart_basis)
self.ra_latt_vec_cart_basis = numpy.array(numpy.reshape(ra_vec, (3, 1)))
self.rb_latt_vec_cart_basis = numpy.array(numpy.reshape(rb_vec, (3, 1)))
self.rc_latt_vec_cart_basis = numpy.array(numpy.reshape(rc_vec, (3, 1)))
[docs] def buildCrystalUnitCell(self, mmcrystal_handle):
"""
Build a crystal unit cell.
:type mmcrystal_handle: handle
:param mmcrystal_handle: the mmcrystal handle
"""
self.crystal_unit_cell = self.asymmetric_unit.copy()
asu_atoms = analyze.evaluate_asl(self.crystal_unit_cell,
'atom.%s' % self.ASU_LABEL_KEY)
# make the unit cell
# note that the following function deletes all properties and
# overwrites the space group as P 1, presumably because once the
# cell is made it loses all symmetries but translational symmetries
# (I choose to keep the space group record as it was originally set
# because everyone knows this fact already)
mm.mmcrystal_generate_cell_ct(mmcrystal_handle, self.crystal_unit_cell)
# overwrite the P 1 space group back to the original space group, so
# that the users can work with different views from the GUI
self.crystal_unit_cell.property[self.SPACE_GROUP_KEY] = self.space_group
# Update crystal Z, needed for the PDB export (SHARED-7022)
self.crystal_unit_cell.property[mm.M2IO_PDB_CRYSTAL_Z] = \
self.asymmetric_unit.property.get(mm.M2IO_PDB_CRYSTAL_Z, 1)
# update coarse grain property
if msutils.is_coarse_grain(self.asymmetric_unit):
coarsegrain.set_as_coarse_grain(self.crystal_unit_cell)
# Preserve original ASU atom property (MATSCI-6914)
for atom in self.crystal_unit_cell.atom:
atom.property[self.ASU_LABEL_KEY] = atom.index in asu_atoms
[docs] def labelSymEquivPos(self):
"""
Label the symmetry equivalent positions.
"""
asu_dim = self.asymmetric_unit.atom_total
# the unit cell is generated by transforms on the original ASU, each
# is called a mate, just unpack them in the correct order
for sym_label in range(1, asu_dim + 1):
for index in range(sym_label, self.crystal_unit_cell.atom_total + 1, \
asu_dim):
self.crystal_unit_cell.atom[index].property[self.SYMMETRY_LABEL_KEY] = \
sym_label
[docs] def labelAsuAtoms(self, astructure):
"""
Label the atoms that make up an ASU, i.e. the symmetry unique atoms.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure whose atoms will be labeled
"""
# Prioritize already present ASU labels, if present (MATSCI-6914)
if analyze.evaluate_asl(astructure, 'atom.%s' % self.ASU_LABEL_KEY):
# Remove SYMMETRY_LABEL_KEY atom property, since these are not used
# anywhere (MATSCI-6914)
msutils.remove_atom_property(astructure, self.SYMMETRY_LABEL_KEY)
return
# loop over atoms and mark the first atom of a given symmetry as an ASU
# atom and the others as non-ASU atoms
sym_labeled = []
for aatom in astructure.atom:
sym_label = aatom.property[self.SYMMETRY_LABEL_KEY]
if sym_label not in sym_labeled:
aatom.property[self.ASU_LABEL_KEY] = True
sym_labeled.append(sym_label)
else:
aatom.property[self.ASU_LABEL_KEY] = False
# Remove SYMMETRY_LABEL_KEY atom property, since these are not used
# anywhere (MATSCI-6914)
msutils.remove_atom_property(astructure, self.SYMMETRY_LABEL_KEY)
[docs] def doPropertyEvaluation(self):
"""
Compute some properties of the unit cell. In order for the
properties to be consistent with their standard definitions the
provided unit cell must be void of the redundant edge, meaning that
the fractionals must be defined on n < f < n + 1 rather than
n < f <= n + 1, as well as any other redundancies.
"""
set_physical_properties(self.crystal_unit_cell)
self.unit_cell_formula = self.crystal_unit_cell.property[
self.UNIT_CELL_FORMULA_KEY]
self.unit_cell_volume = self.crystal_unit_cell.property[
self.UNIT_CELL_VOLUME_KEY]
self.unit_cell_density = self.crystal_unit_cell.property[
self.UNIT_CELL_DENSITY_KEY]
[docs] def setStructureProperties(self):
"""
Set some structure properties.
"""
# some third-party software may require the space group ID so make
# that a structure property as well
self.crystal_unit_cell.property[self.SPACE_GROUP_ID_KEY] = \
self.spgobj.space_group_id
# honor incoming settings, for example like 'anchor_%s_%s_%s' % (0, 0, 0)
# to fix the origin, etc. but use 'center_structure' if absent
pbc_position = self.asymmetric_unit.property.get(
PBC_POSITION_KEY, CENTER_PBC_POSITION)
self.crystal_unit_cell.property[PBC_POSITION_KEY] = pbc_position
[docs] def setChorusProperties(self, astructure, ncella=1, ncellb=1, ncellc=1):
"""
Set the chorus structure properties on the given structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure that needs the chorus
properties
:type ncella: int
:param ncella: the number of cells along a
:type ncellb: int
:param ncellb: the number of cells along b
:type ncellc: int
:param ncellc: the number of cells along c
"""
lengths = [
ncella * self.a_param, ncellb * self.b_param, ncellc * self.c_param
]
params = lengths + [self.alpha_param, self.beta_param, self.gamma_param]
for key, value in zip(self.CHORUS_BOX_KEYS,
get_chorus_from_params(params)):
astructure.property[key] = value
[docs] def printCrystalParams(self):
"""
Print some crystal parameters.
"""
def format_vector(vector):
new_vector = []
for coord in vector.flatten():
element = '%.4f' % round(coord, self.NUM_DECIMAL_COORDS)
new_vector.append(element)
new_vector = ', '.join(new_vector)
new_vector = '[' + new_vector + ']'
return new_vector
HEADER = 'Crystal parameters'
unit_cell_formula = textlogger.get_param_string('Unit cell formula',
self.unit_cell_formula,
self.MSGWIDTH)
unit_cell_volume = '%.4f' % round(self.unit_cell_volume,
self.NUM_DECIMAL_COORDS)
unit_cell_volume = textlogger.get_param_string(
'Unit cell volume / Ang.^3', unit_cell_volume, self.MSGWIDTH)
unit_cell_density = '%.4f' % round(self.unit_cell_density,
self.NUM_DECIMAL_COORDS)
unit_cell_density = textlogger.get_param_string(
'Unit cell density / g/cm^3', unit_cell_density, self.MSGWIDTH)
definition_id = textlogger.get_param_string('Definition ID',
self.spgobj.definition_id,
self.MSGWIDTH)
space_group_id = textlogger.get_param_string('Space group ID',
self.spgobj.space_group_id,
self.MSGWIDTH)
ichoice = textlogger.get_param_string('Index of setting choice',
self.spgobj.ichoice,
self.MSGWIDTH)
num_choices = textlogger.get_param_string('Number of choices',
self.spgobj.num_choices,
self.MSGWIDTH)
spg_setting = textlogger.get_param_string('Space Group Setting',
self.spgobj.spg_setting,
self.MSGWIDTH)
short_name = textlogger.get_param_string( \
'Space group short Hermann-Mauguin symbol',
self.spgobj.space_group_short_name, self.MSGWIDTH)
full_name = textlogger.get_param_string( \
'Space group full Hermann-Mauguin symbol',
self.spgobj.space_group_full_name, self.MSGWIDTH)
point_group = textlogger.get_param_string('Point group',
self.spgobj.point_group_name,
self.MSGWIDTH)
ncent = textlogger.get_param_string('Number of centering operations',
self.spgobj.num_centering_opers,
self.MSGWIDTH)
nprim = textlogger.get_param_string('Number of primary operations',
self.spgobj.num_primary_opers,
self.MSGWIDTH)
nsym = textlogger.get_param_string('Number of symmetry operations',
self.spgobj.num_symmetry_opers,
self.MSGWIDTH)
crystal_system = textlogger.get_param_string(
'Crystal system', self.spgobj.crystal_system.name, self.MSGWIDTH)
length_a = '%.4f' % round(self.a_param, self.NUM_DECIMAL_COORDS)
length_a = textlogger.get_param_string('Length lattice vector a / Ang.',
length_a, self.MSGWIDTH)
length_b = '%.4f' % round(self.b_param, self.NUM_DECIMAL_COORDS)
length_b = textlogger.get_param_string('Length lattice vector b / Ang.',
length_b, self.MSGWIDTH)
length_c = '%.4f' % round(self.c_param, self.NUM_DECIMAL_COORDS)
length_c = textlogger.get_param_string('Length lattice vector c / Ang.',
length_c, self.MSGWIDTH)
angle_alpha = '%.4f' % round(self.alpha_param, self.NUM_DECIMAL_COORDS)
angle_alpha = textlogger.get_param_string(
'Lattice angle alpha / Degree', angle_alpha, self.MSGWIDTH)
angle_beta = '%.4f' % round(self.beta_param, self.NUM_DECIMAL_COORDS)
angle_beta = textlogger.get_param_string('Lattice angle beta / Degree',
angle_beta, self.MSGWIDTH)
angle_gamma = '%.4f' % round(self.gamma_param, self.NUM_DECIMAL_COORDS)
angle_gamma = textlogger.get_param_string(
'Lattice angle gamma / Degree', angle_gamma, self.MSGWIDTH)
ncell_a = textlogger.get_param_string('Number of unit cells along a',
self.ncella, self.MSGWIDTH)
ncell_b = textlogger.get_param_string('Number of unit cells along b',
self.ncellb, self.MSGWIDTH)
ncell_c = textlogger.get_param_string('Number of unit cells along c',
self.ncellc, self.MSGWIDTH)
origin = textlogger.get_param_string('Origin', self.origin,
self.MSGWIDTH)
a_latt_vec_cart_basis = textlogger.get_param_string( \
'Lattice vector a in Cartesian basis',
format_vector(self.a_latt_vec_cart_basis), self.MSGWIDTH)
b_latt_vec_cart_basis = textlogger.get_param_string( \
'Lattice vector b in Cartesian basis',
format_vector(self.b_latt_vec_cart_basis), self.MSGWIDTH)
c_latt_vec_cart_basis = textlogger.get_param_string( \
'Lattice vector c in Cartesian basis',
format_vector(self.c_latt_vec_cart_basis), self.MSGWIDTH)
do_bonding = textlogger.get_param_string('Do bonding', self.do_bonding,
self.MSGWIDTH)
do_bond_orders = textlogger.get_param_string('Do bond orders',
self.do_bond_orders,
self.MSGWIDTH)
do_translation = textlogger.get_param_string('Do translation',
self.do_translation,
self.MSGWIDTH)
num_sym_unique = textlogger.get_param_string(
'Number of symmetry unique atoms', self.num_sym_unique,
self.MSGWIDTH)
self.logger.info(HEADER)
self.logger.info('-' * len(HEADER))
self.logger.info('')
self.logger.info(unit_cell_formula)
self.logger.info(unit_cell_volume)
self.logger.info(unit_cell_density)
self.logger.info('')
self.logger.info(definition_id)
self.logger.info(space_group_id)
self.logger.info(ichoice)
self.logger.info(num_choices)
self.logger.info(spg_setting)
self.logger.info(short_name)
self.logger.info(full_name)
self.logger.info(point_group)
self.logger.info(ncent)
self.logger.info(nprim)
self.logger.info(nsym)
self.logger.info(crystal_system)
self.logger.info('')
self.logger.info(length_a)
self.logger.info(length_b)
self.logger.info(length_c)
self.logger.info(angle_alpha)
self.logger.info(angle_beta)
self.logger.info(angle_gamma)
self.logger.info(ncell_a)
self.logger.info(ncell_b)
self.logger.info(ncell_c)
self.logger.info('')
self.logger.info(origin)
self.logger.info(a_latt_vec_cart_basis)
self.logger.info(b_latt_vec_cart_basis)
self.logger.info(c_latt_vec_cart_basis)
self.logger.info('')
self.logger.info(do_bonding)
self.logger.info(do_bond_orders)
self.logger.info(do_translation)
self.logger.info(num_sym_unique)
self.logger.info('')
[docs] def setAsuAtomsFalse(self, astructure):
"""
Set the ASU atom labels to false for this structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object to be updated
"""
for aatom in astructure.atom:
aatom.property[self.ASU_LABEL_KEY] = False
[docs] def buildCrystalSuperCell(self, ncella, ncellb, ncellc):
"""
Build a crystal super cell.
:type ncella: int
:param ncella: the number of unit cells to generate
along lattice vector a
:type ncellb: int
:param ncellb: the number of unit cells to generate
along lattice vector b
:type ncellc: int
:param ncellc: the number of unit cells to generate
along lattice vector c
"""
self.crystal_super_cell = self.crystal_unit_cell.copy()
if ncella == ncellb == ncellc == 1:
return
pbc_bonds = self.pbc_bonds_by_cell[PBCBonds.FIRST_CELL]
blank_cell = self.crystal_super_cell.copy()
self.setAsuAtomsFalse(blank_cell)
unit_cell_dim = self.crystal_super_cell.atom_total
for indexa in range(1, ncella + 1):
for indexb in range(1, ncellb + 1):
for indexc in range(1, ncellc + 1):
cell_indices = (indexa, indexb, indexc)
if indexa == indexb == indexc == 1:
continue
else:
pos = (indexa - 1) * self.a_latt_vec_cart_basis + \
(indexb - 1) * self.b_latt_vec_cart_basis + \
(indexc - 1) * self.c_latt_vec_cart_basis
posx, posy, posz = pos
next_cell = blank_cell.copy()
transform.translate_structure(next_cell, posx, posy,
posz)
# save atomic indices of new cell, if redundant boundaries
# are wanted then we must delete any overlapping boundary
# atoms of the new cell
dim_so_far = self.crystal_super_cell.atom_total
indicies = list(
range(dim_so_far + 1,
dim_so_far + unit_cell_dim + 1))
self.crystal_super_cell.extend(next_cell)
num_deleted_atoms = 0
next_pbc_bonds = pbc_bonds.getOffsetPBCBonds(cell_indices, \
dim_so_far - num_deleted_atoms)
self.pbc_bonds_by_cell[cell_indices] = next_pbc_bonds
self.indicies_by_cell[cell_indices] = list(indicies)
# due to the way that we are growing the cell it is more
# efficient to bond ahead only for the last cells
if self.do_bonding:
next_pbc_bonds.connectToCells(
self.crystal_super_cell, ncella, ncellb, ncellc,
PBCBond.BEHIND)
if indexa == ncella or indexb == ncellb or indexc == ncellc:
next_pbc_bonds.connectToCells(
self.crystal_super_cell, ncella, ncellb,
ncellc, PBCBond.AHEAD)
[docs] def setOrigin(self):
"""
Set the origin.
"""
if self.origin is None:
pbc_position = self.asymmetric_unit.property.get(PBC_POSITION_KEY)
if pbc_position and pbc_position.startswith(ANCHOR_PREFIX):
xyzs = get_carts_from_anchor_string(pbc_position)
fracts = numpy.dot(self.xyz2cry, numpy.array(xyzs))
self.origin = fracts.tolist()
else:
self.origin = ParserWrapper.ORIGIN
[docs] def orchestrate(self):
"""
Orchestrate the construction of the crystal.
"""
# do some preliminary set up
self.getLatticeParameters()
self.checkInputParams()
self.updateLatticeProperties()
# get the mmcrystal handle and fractional to-from Cartesian
# transforms
mmcrystal_handle = mm.mmcrystal_new()
try:
self._setCrystalSymmetryFromHandle(mmcrystal_handle)
self.determineBasisVectors()
# set the origin
self.setOrigin()
# check redundancies in fractional coordinates, see MATSCI-4633
# where it was decided to perform this check only for non-P1 or
# non-Desmond cells
state = self.asymmetric_unit.property[self.SPACE_GROUP_KEY] == \
P1_SPACE_GROUP_SYMBOL or self.asymmetric_unit.property.get(
CT_TYPE) == CT_TYPE.VAL.FULL_SYSTEM
if not state:
asu_copy = self.asymmetric_unit.copy()
CheckInput().checkFractionalRedundancies(asu_copy, self.logger)
self.handleSettings()
# build the unit cell and label the symmetry equivalent atoms
self.buildCrystalUnitCell(mmcrystal_handle)
finally:
mm.mmcrystal_delete(mmcrystal_handle)
self.labelSymEquivPos()
# The function mmcrystal_generate_cell_ct which applies the symmetry
# operations of the space group to the ASU to provide the unit cell
# treats all positions as general positions. This means that if the
# ASU features a special position that there will be duplicate copies
# of the atom at that position which are exactly coincident with each
# other. These duplicates need to be removed. Note that the mmcrystal
# code contains a function called mmcrystal_remove_specialpos which is
# supposed to do this however it works only within the mates option,
# not the CT option, is buggy and crashes upon application, has
# never been called, and appears to be difficult to SWIG wrap, also
# we need the unit cell void of all redundancies so that we obtain
# valid properties
# Preserve bonding only if bonding/bond orders recalculation wasn't
# requested (it can be slow MATSCI-7196)
preserve_bonding = not (self.do_bonding or self.do_bond_orders)
if self.overlap_tresh:
delete_duplicate_atoms(self.crystal_unit_cell,
fract_offset=self.fract_offset,
duplicate_thresh=self.overlap_tresh,
preserve_bonding=preserve_bonding)
self.setChorusProperties(self.crystal_unit_cell)
self.doPropertyEvaluation()
self.setStructureProperties()
# now that we are finished preparing the unit cell label those atoms that
# define an ASU, i.e. symmetry unique atoms, and set the number of symmetry
# unique atoms, this is not simply the number of atoms in the input ASU
# as it may have included redundant atoms
self.labelAsuAtoms(self.crystal_unit_cell)
self.num_sym_unique = len([aatom.index for aatom in self.crystal_unit_cell.atom \
if aatom.property.get(self.ASU_LABEL_KEY)])
# do bonding, since it is the unit cell we want to use the unit cell
# chorus properties, always create PBC bonds because we will use them
# for the crystal build, depending on self.do_bonding create regular
# bonds as well (it is possible that a user might want PBC bonds but
# not the regular bonds to be assigned as the later might already be
# there, for example as in a protein, etc.), possibly remove the PBC
# bonds later
if not self.use_existing_pbc_bonds:
if self.do_bonding:
maximally_bonded_atoms, bonds_dict, pbc_bonds_dict = \
connect_atoms(self.crystal_unit_cell, cov_offset=self.cov_offset,
delete_existing=True, pbc_bonding=True, only_pbc_bonds=False,
max_valencies=self.max_valencies)
else:
maximally_bonded_atoms, bonds_dict, pbc_bonds_dict = [], {}, {}
if maximally_bonded_atoms and self.logger:
warn_msg = ('The following atoms have a maximal number of '
'bonds and thus cannot be further bonded: %s.') % \
maximally_bonded_atoms
self.logger.warning(warn_msg)
else:
pbc_bonds_dict = get_pbc_bonds_dict(self.crystal_unit_cell)
tmp_st = self.crystal_unit_cell.copy()
clusterstruct.contract_structure(tmp_st)
label_pbc_bonds(tmp_st)
tmp_pbc_bonds_dict = get_pbc_bonds_dict(tmp_st)
if not pbc_bonds_dict or pbc_bonds_dict and not tmp_pbc_bonds_dict:
self.crystal_unit_cell = tmp_st
pbc_bonds_dict = tmp_pbc_bonds_dict
pbc_bonds = PBCBonds(pbc_bonds_dict=pbc_bonds_dict,
cell_size=self.crystal_unit_cell.atom_total,
cov_offset=self.cov_offset)
# set chorus properties for the super cell
self.setChorusProperties(self.crystal_unit_cell, \
ncella=self.ncella, ncellb=self.ncellb, ncellc=self.ncellc)
# do bond orders
if self.do_bond_orders:
self.crystal_unit_cell = assign_bond_orders(self.crystal_unit_cell,
logger=self.logger)
pbc_bonds.updatePBCBondOrders(self.crystal_unit_cell)
# set the cell deltas for the PBC bonds
lattice_vectors = [self.a_latt_vec_cart_basis, \
self.b_latt_vec_cart_basis, self.c_latt_vec_cart_basis]
pbc_bonds.setDeltasToNeighboringCells(self.crystal_unit_cell, \
lattice_vectors)
# finally print some info
if self.logger:
self.printCrystalParams()
# update collections
self.pbc_bonds_by_cell[(1, 1, 1)] = pbc_bonds
self.indicies_by_cell[(1, 1, 1)] = \
list(range(1, self.crystal_unit_cell.atom_total + 1))
# build the super cell via translations of the unit cell
self.buildCrystalSuperCell(self.ncella, self.ncellb, self.ncellc)
# if the unit cell contains a diatomic molecule that has PBC bonds
# and we are doing extents and assigning bond orders then we will
# need to assign bond orders to the entire super cell to catch
# resonances, see MATSCI-2032
extents = [self.ncella, self.ncellb, self.ncellc]
extents_present = any((x > 1 for x in extents))
if self.do_bond_orders and extents_present:
for amol in self.crystal_unit_cell.molecule:
atoms = amol.getAtomIndices()
if len(atoms) == 2:
abond = self.crystal_unit_cell.getBond(*atoms)
if abond.property.get(PBCBond.PBC_BOND_KEY):
self.crystal_super_cell = \
assign_bond_orders(self.crystal_super_cell, logger=self.logger)
break
# if a translation back to the first cell was called for then do it
if self.do_translation:
translate_to_cell(self.crystal_super_cell,
fract_offset=self.fract_offset,
origin=self.origin,
extents=extents)
elif self.translate_centroids:
translate_molecules(self.crystal_super_cell,
centroids=True,
fract_offset=self.fract_offset)
if self.do_translation or self.translate_centroids:
label_pbc_bonds(self.crystal_super_cell)
[docs]def get_reciprocal_lattice_vectors(a_vec, b_vec, c_vec):
"""
Return the reciprocal lattice vectors.
: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: numpy.array
:return: the three reciprocal lattice vectors
"""
def get_reciprocal_vec(vec1, vec2, vec3):
vec23 = numpy.cross(vec2, vec3)
return old_div(vec23, numpy.dot(vec1, vec23))
ra_vec = get_reciprocal_vec(a_vec, b_vec, c_vec)
rb_vec = get_reciprocal_vec(b_vec, c_vec, a_vec)
rc_vec = get_reciprocal_vec(c_vec, a_vec, b_vec)
return numpy.array([ra_vec, rb_vec, rc_vec])
[docs]def get_collapsed_index(abc, alimit, blimit, climit):
"""
Given a three dimensional grid of integers defined on
[1, limit] for the given a, b, and c limits and a traversal
path of c then b then a return the number of integers
traversed in order to reach the given abc integer index
triple.
:type abc: tuple
:param abc: the integer indices
:type alimit: int
:param alimit: the upper bound on a (unused)
:type blimit: int
:param blimit: the upper bound on b
:type climit: int
:param climit: the upper bound on c
:rtype: int
:return: the number of integers traversed
"""
a_index, b_index, c_index = abc
return (a_index - 1) * blimit * climit + (b_index - 1) * climit + (c_index -
1) + 1
[docs]def modified_sawtooth(n, x):
"""
Given a positive integer variable x in [0, n+1] return a signal
from a modified sawtooth function. This function is linear on
[1, n] but is n for x = 0 and is 1 for x = n+1.
:type n: int
:param n: the period
:type x: int
:param x: the variable
:rtype: int
:return: the signal
"""
if x == 0:
return n
elif x == n + 1:
return 1
else:
return x
[docs]def assign_bond_orders(astructure, logger=None):
"""
Return a copy of the input structure that has bond
orders assigned.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object for which bond
orders will be assigned
:type logger: logging.getLogger or None
:param logger: output logger or None if there isn't one
:rtype: schrodinger.Structure.structure
:return: a copy of the input structure with the bond
orders assigned
"""
# do this by passing nothing other than individual
# molecule structures, extracted from the given
# structure, to assignbondorders.assign_st, as this
# is more efficient than passing the entire structure
# and either assigning all atoms or atom lists for
# molecules
bo_astructure = structure.create_new_structure()
bo_astructure.property = dict(astructure.property)
reorder_map = []
warn = False
for amol in astructure.molecule:
amol_old_indices = sorted(amol.getAtomIndices())
amol_natoms = len(amol_old_indices)
offset = bo_astructure.atom_total + 1
amol_new_indices = list(range(offset, amol_natoms + offset))
reorder_map.extend(list(zip(amol_old_indices, amol_new_indices)))
amol_astructure = amol.extractStructure(copy_props=True)
with ioredirect.IOCapture() as obj:
assignbondorders.assign_st(amol_astructure)
msg = obj.getvalue().strip()
if msg:
warn = True
bo_astructure.extend(amol_astructure)
if warn and logger:
msg = ('The bond order assignment may have been unsuccessful '
'for at least one of the molecules in the system.')
logger.warning(msg)
old_idxs, new_idxs = list(zip(*reorder_map))
reorder_list = [old_idxs.index(new_idx) + 1 for new_idx in new_idxs]
return build.reorder_atoms(bo_astructure, reorder_list)
[docs]def assign_bond_orders_w_mmlewis(astructure, fix_metals=True, logger=None):
"""
Return a copy of the input structure that has bond
orders assigned.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure object for which bond
orders will be assigned
:type fix_metals: bool
:param fix_metals: fix metals coming from mmlewis
:type logger: logging.getLogger or None
:param logger: output logger or None if there isn't one
:rtype: schrodinger.Structure.structure
:return: a copy of the input structure with the bond
orders assigned
"""
mm.mmlewis_initialize(mm.MMERR_DEFAULT_HANDLER)
old_handler = mm.mmlewis_get_error_handler()
new_handler = ErrorHandler(queued=True, silent=True)
new_handler.push_level(6)
mm.mmlewis_set_error_handler(new_handler.handle)
# do this by passing nothing other than individual
# molecule structures, extracted from the given
# structure, to mmlewis_apply, as this is more
# efficient than passing the entire structure
warn = False
bo_astructure = structure.create_new_structure()
bo_astructure.property = dict(astructure.property)
reorder_map = []
for amol in astructure.molecule:
amol_old_indices = sorted(amol.getAtomIndices())
amol_natoms = len(amol_old_indices)
offset = bo_astructure.atom_total + 1
amol_new_indices = list(range(offset, amol_natoms + offset))
reorder_map.extend(list(zip(amol_old_indices, amol_new_indices)))
amol_astructure = amol.extractStructure(copy_props=True)
try:
mm.mmlewis_apply(amol_astructure)
except:
warn = True
bo_astructure.extend(amol_astructure)
old_idxs, new_idxs = list(zip(*reorder_map))
reorder_list = [old_idxs.index(new_idx) + 1 for new_idx in new_idxs]
bo_astructure = build.reorder_atoms(bo_astructure, reorder_list)
if warn and logger:
msg = ('The bond order assignment was unsuccessful for at '
'least one of the molecules in the system.')
logger.warning(msg)
mm.mmlewis_set_error_handler(old_handler)
mm.mmlewis_terminate()
if fix_metals:
fix_metal_bonding(bo_astructure)
return bo_astructure
[docs]def get_chorus_properties(astructure):
"""
Return a tuple containing the nine chorus properties
of the given structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure that has the chorus
properties defined
:rtype: tuple
:return: contains the nine chorus properties
"""
return tuple((astructure.property[x] for x in Crystal.CHORUS_BOX_KEYS))
[docs]def get_params_from_chorus(chorus_properties):
"""
Return the a, b, c, alpha, beta, and gamma lattice properties
from the nine chorus properties.
:type chorus_properties: list
:param chorus_properties: contains the nine chorus properties,
i.e. ax, ay, az, bx, ..., cz
:rtype: list
:return: a, b, c, alpha, beta, and gamma lattice properties
"""
vectors = [numpy.array(chorus_properties[x:x + 3]) for x in [0, 3, 6]]
return get_params_from_vectors(*vectors)
[docs]def get_chorus_from_params(params):
"""
Return the nine chorus properties, i.e. [ax, ay, az, bx, ..., cz],
from the six lattice parameters a, b, c, alpha, beta, and gamma.
:type params: list
:param params: contains the a, b, c, alpha, beta, and gamma lattice
parameters
:rtype: list
:return: contains the nine chorus properties
"""
a_vec, b_vec, c_vec = get_lattice_vectors(*params)
return list(a_vec) + list(b_vec) + list(c_vec)
[docs]def get_volume_from_params(params):
"""
Return cell volume (in Angstrom^3) from cell parameters.
:type params: list
:param params: contains the a, b, c, alpha, beta, and gamma lattice
parameters
:rtype: float
:return: Cell volume in Angstrom^3
"""
a_vec, b_vec, c_vec = get_lattice_vectors(*params)
vol = get_volume_from_vecs([a_vec, b_vec, c_vec])
return vol
[docs]def get_volume_from_vecs(vecs):
"""
Return cell volume (in Angstrom^3) from lattice vectors.
:type params: list
:param params: lattice vectors
parameters
:rtype: float
:return: Cell volume in Angstrom^3
"""
return abs(numpy.dot(numpy.cross(vecs[0], vecs[1]), vecs[2]))
[docs]def get_cell_pairs(astructure,
cell_distance,
pbc_bonding=True,
atom_indices=None,
chorus_properties=None):
"""
Using a distance cell that optionally honors a PBC return
a list of tuples of atom index pairs that are within the specified
distance.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure containing the pairs
:type cell_distance: float
:param cell_distance: the distance used in the distance cell,
if using a PBC then min([cell_distance, a, b, c]) is actually
what is used
:type pbc_bonding: bool
:param pbc_bonding: if True and the chorus box properties exist
on the incoming structure then the distance cell will honor
the PBC, otherwise any PBC is not considered
:type atom_indices: list or None
:param atom_indices: a list of atom indices to search for pairs,
each atom is searched for pairs that may or may not already be in
this list, if None then all atoms are searched
:type chorus_properties: list or None
:param chorus_properties: contains the nine chorus box properties
that define a PBC, i.e. [ax, ay, az, bx, ..., cz], if None then
the chorus structure properties will be used if available and if
not available then no PBC will be used in the DistanceCell
:rtype: list
:return: contains tuples of unique atom index pairs within the
specified distance
"""
if not chorus_properties:
try:
chorus_properties = get_chorus_properties(astructure)
except KeyError:
chorus_properties = None
if chorus_properties and pbc_bonding:
pbc = infrastructure.PBC(*chorus_properties)
cell, _ = clusterstruct.create_distance_cell(astructure,
cell_distance,
pbc=pbc)
else:
cell = infrastructure.DistanceCell(astructure, cell_distance)
if not atom_indices:
atom_indices = list(range(1, astructure.atom_total + 1))
# Dict that hold pairs as keys and minimum distances as values
all_pairs = {}
for a_index in atom_indices:
a_atom = astructure.atom[a_index]
for match in cell.query_atoms(*(a_atom.xyz)):
if match.getIndex() == a_index:
continue
pair = frozenset([a_index, match.getIndex()])
dist = math.sqrt(match.getDistanceSquared())
# After SHARED-7697, more than one atom image can be present in
# matches. Keep only shortest distance
if pair not in all_pairs or all_pairs[pair] > dist:
all_pairs[pair] = dist
# Sort by distance, return only keys (pairs)
all_pairs = [
tuple(sorted(k))
for k, _ in sorted(all_pairs.items(), key=lambda x: x[1])
]
return all_pairs
[docs]class PBCBond(object):
"""
Class to manage a PBC bond, i.e. a long bond connecting
two real atoms on opposite sides of a unit cell that is
used in lieu of the bond between the real and image atoms.
"""
PBC_BOND_THRESHOLD = 1e-3
PBC_BOND_KEY = mm.M2IO_DATA_PBC_BOND
ALSO_REG_BOND_KEY = mm.M2IO_DATA_PBC_ALSO_REG_BOND
PBC_BOND_COLOR_KEY = mm.M2IO_DATA_COLOR
# this is blue
PBC_BOND_COLOR = 4
AHEAD = 'ahead'
BEHIND = 'behind'
[docs] def __init__(self, atom1, atom2, order, also_reg_bond):
"""
Create an instance.
:type atom1: int
:param atom1: the first atom index of the bond
:type atom2: int
:param atom2: the second atom index of the bond
:type order: int
:param order: the bond order
:type also_reg_bond: bool
:param also_reg_bond: indicates whether this PBC bond
is also a regular bond, meaning that one of the atoms
of the PBC bond is covalently bound to two copies of
the other atom, one inside the cell and one outside the
cell
"""
self.atom1 = atom1
self.atom2 = atom2
self.order = order
self.also_reg_bond = also_reg_bond
self.cell_deltas_ahead = {}
self.cell_deltas_behind = {}
def __repr__(self):
"""
Define a class representation.
:rtype: str
:return: string representation of the class
"""
msg = '(%s, %s): %s\nAlso regular bond: %s\nAhead: %s\nBehind: %s\n\n'
return msg % (self.atom1, self.atom2, self.order, self.also_reg_bond, \
self.cell_deltas_ahead, self.cell_deltas_behind)
[docs] def setDeltasToNeighboringCells(self, atom1_vec, atom2_vec, \
spanning_vectors):
"""
Set two dictionaries, one for moving to neighboring cells
ahead and one for moving behind. Keys are tuples of
integer cell deltas and values are (tail, head) tuples
giving a directionality of this PBC bond along the given
direction.
:type atom1_vec: numpy.array
:param atom1_vec: vector to the first atom of the PBC bond
:type atom2_vec: numpy.array
:param atom2_vec: vector to the second atom of the PBC bond
:type spanning_vectors: dict
:param spanning_vectors: keys are tuples of cell index triples,
values are tuples of normalized spanning vectors for the cell
and their original lengths
"""
pbc_bond_vec = atom2_vec - atom1_vec
cell_deltas = []
for triple, (vec, alen) in spanning_vectors.items():
projection = numpy.dot(pbc_bond_vec, vec)
if abs(projection) < 0.5 * alen - self.PBC_BOND_THRESHOLD:
continue
anti = projection < 0
cell_deltas.append((triple, anti))
for (cell_delta, anti) in cell_deltas:
rev_cell_delta = tuple([-1 * idx for idx in cell_delta])
pair = (self.atom1, self.atom2)
if anti:
pair = pair[::-1]
self.cell_deltas_ahead[cell_delta] = pair
self.cell_deltas_behind[rev_cell_delta] = pair[::-1]
[docs] def getNeighborPBCBond(self, cell_indices, ncella, ncellb, ncellc,
cell_size, cell_delta, tail_head):
"""
Return a (tail, head) ordered tuple of atom indices for
the neighboring PBC bond in the cell given by the cell delta.
If this cell is the first or last cell then it will wrap
around to the ending or beginning cell, respectively.
:type cell_indices: tuple
:param cell_indices: a triple of cell indices
:type ncella: int
:param ncella: the number of cells along a
:type ncellb: int
:param ncellb: the number of cells along b
:type ncellc: int
:param ncellc: the number of cells along c
:type cell_size: int
:param cell_size: the number of atoms in a cell
:type cell_delta: tuple
:param cell_delta: a triple of cell deltas which provide
the neighboring cells location
:type tail_head: tuple
:param tail_head: provides the tail and head atom indices
for this PBC bond given the cell delta
:rtype: tuple
:return: atom indices for the neighboring PBC bond
"""
extents = [ncella, ncellb, ncellc]
next_cell = tuple(map(sum, list(zip(cell_indices, cell_delta))))
next_cell = tuple(
(modified_sawtooth(*x) for x in zip(extents, next_cell)))
# natom_diff is signed
natom_diff = \
get_natom_btw_two_cells(next_cell, cell_indices, \
extents, cell_size)
return tuple((x - natom_diff for x in tail_head))
[docs] def offsetAtomData(self, offset):
"""
Offset the atom data of this instance.
:type offset: int
:param offset: the offset to use
"""
def update_dict(adict):
for key in list(adict):
value = adict[key]
value = tuple((x + offset for x in value))
adict[key] = value
return adict
self.atom1 += offset
self.atom2 += offset
self.cell_deltas_ahead = update_dict(self.cell_deltas_ahead)
self.cell_deltas_behind = update_dict(self.cell_deltas_behind)
[docs]class PBCBonds(object):
"""
Class to manage a collection of PBC bonds.
"""
FIRST_CELL = (1, 1, 1)
[docs] def __init__(self,
cell_indices=None,
pbc_bonds_dict=None,
cell_size=None,
cov_offset=ParserWrapper.COV_OFFSET):
"""
Create an instance.
:type cell_indices: tuple or None
:param cell_indices: a triple of cell indices indicating
to which cell the given PBC bonds belong, if None then the
first cell will be used
:type pbc_bonds_dict: dict or None
:param pbc_bonds_dict: used to create a dictionary of PBCBond
objects, keys are tuples of PBC bond atom index pairs,
values are tuples of bond orders and booleans indicating whether
the PBC bond is also a regular bond, if None then the dictionary
of PBCBond ojects will be empty
:type cell_size: int
:param cell_size: the number of atoms in a cell
:type cov_offset: float
:param cov_offset: the maximum distance for a connection is
the sum of the covalent radii of the two atoms weighted by
cov_factor plus this offset in angstrom, increasing this
value will increase the number of connections
"""
if cell_indices is None:
cell_indices = self.FIRST_CELL
self.cell_indices = cell_indices
self.pbc_bonds = {}
if pbc_bonds_dict:
self.setPBCBonds(pbc_bonds_dict)
self.cell_size = cell_size
self.cov_offset = cov_offset
[docs] def setPBCBonds(self, pbc_bonds_dict):
"""
Create a dictionary of PBCBond objects from a dictionary
containing PBC bonds.
:type pbc_bonds_dict: dict
:param pbc_bonds_dict: keys are tuples of PBC bond atom
index pairs, values are tuples of bond orders and booleans
indicating whether the PBC bond is also a regular bond
"""
for pair, (order, also_reg_bond) in pbc_bonds_dict.items():
self.pbc_bonds[pair] = \
PBCBond(pair[0], pair[1], order, also_reg_bond)
[docs] def updatePBCBondOrders(self, astructure):
"""
Update the bond orders in this instance given a
structure with updated PBC bond orders.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure with the updated PBC
bond orders
"""
if not self.pbc_bonds:
return
for abond in astructure.bond:
key = tuple(sorted([abond.atom1.index, abond.atom2.index]))
if abond.property.get(PBCBond.PBC_BOND_KEY):
self.pbc_bonds[key].order = abond.order
def __repr__(self):
"""
Define a class representation.
:rtype: str
:return: string representation of the class
"""
msg = 'Cell: ' + str(self.cell_indices) + '\n\n'
for abond in self.pbc_bonds.values():
msg += repr(abond)
return msg
[docs] def setDeltasToNeighboringCells(self, astructure, lattice_vectors):
"""
For each PBC bond in this cell determine the cell deltas
that are needed to move ahead and behind to relevant
neighboring cells.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure with the PBC bonds
:type lattice_vectors: list of numpy.array
:param lattice_vectors: the lattice a, b, and c vectors
"""
spanning_vectors = get_normalized_spanning_vectors(lattice_vectors)
for pair, pbc_bond in self.pbc_bonds.items():
atom1_vec, atom2_vec = \
[numpy.array(astructure.atom[idx].xyz) for idx in pair]
pbc_bond.setDeltasToNeighboringCells(atom1_vec, atom2_vec, \
spanning_vectors)
[docs] def cleanUpPBCBonds(self, astructure, pairs, delete):
"""
Clean up the specified PBC bonds in the structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure with the PBC bonds
:type pairs: list
:param pairs: tuples of PBC bonding atom index pairs
:type delete: bool
:param delete: if True then the PBC bonds will be deleted
"""
for pair in pairs:
abond = astructure.getBond(*pair)
if abond:
abond.property.pop(PBCBond.PBC_BOND_KEY, None)
abond.property.pop(PBCBond.ALSO_REG_BOND_KEY, None)
if delete:
abond.delete()
[docs] def getOffsetPBCBonds(self, cell_indices, offset):
"""
Return a PBCBonds instance for the provided cell indices in
which the atom indices have been offset by the given amount.
:type cell_indices: tuple
:param cell_indices: a triple of cell indices indicating
to which cell the offset PBC bonds belong
:type offset: int
:param offset: the offset to use in setting the atom indices
:rtype: PBCBonds
:return: a PBCBonds object for the given cell containing
the offset atoms
"""
pbc_bonds_obj = copy.deepcopy(self)
pbc_bonds_obj.cell_indices = cell_indices
pbc_bonds_dict = {}
for pair, pbc_bond in pbc_bonds_obj.pbc_bonds.items():
offset_pair = \
tuple((x + offset for x in pair))
pbc_bond.offsetAtomData(offset)
pbc_bonds_dict[offset_pair] = pbc_bond
pbc_bonds_obj.pbc_bonds = pbc_bonds_dict
return pbc_bonds_obj
[docs] def connectToCells(self, astructure, ncella, ncellb, ncellc, direction):
"""
Connect the PBC bonds in this cell with those relevant
neighboring cells ahead of or behind this one so as to
create real bonds from pairs of PBC bonds or new PBC
bonds. If this cell is the first or last along a given
direction then if possible it will create new PBC bonds
by wrapping around to ending or beginning cells,
respectively.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure for which the connections
are sought
:type ncella: int
:param ncella: the number of cells along a
:type ncellb: int
:param ncellb: the number of cells along b
:type ncellc: int
:param ncellc: the number of cells along c
:type direction: str
:param direction: the direction in which to go for the
next PBC bond, either PBCBond.AHEAD or PBCBond.BEHIND
"""
for pbc_bond in self.pbc_bonds.values():
if direction == PBCBond.AHEAD:
cell_deltas = pbc_bond.cell_deltas_ahead
else:
cell_deltas = pbc_bond.cell_deltas_behind
neighbors_done = set()
for cell_delta, curr_tail_head in cell_deltas.items():
neighbor_tail_head = \
pbc_bond.getNeighborPBCBond(self.cell_indices, \
ncella, ncellb, ncellc, self.cell_size, cell_delta, \
curr_tail_head)
# handle wrapping
neighbor_tail_head_sorted = tuple(sorted(neighbor_tail_head))
if neighbor_tail_head == curr_tail_head or \
neighbor_tail_head_sorted in neighbors_done:
continue
neighbors_done.add(neighbor_tail_head_sorted)
# the index error means that the structure doesn't yet
# have atoms in the necessary cell
try:
pair_objs = [astructure.atom[x] for x in neighbor_tail_head]
except IndexError:
continue
curr_tail, curr_head = curr_tail_head
neighbor_tail, neighbor_head = neighbor_tail_head
pairs = [(curr_head, neighbor_tail), (curr_tail, neighbor_head)]
# Allow maximum number of bonds instead of using default valencies
max_valencies = elementalprops.get_mmct_max_valencies()
bonds_made = False
for pair in pairs:
maximally_bonded_atoms, bonds_dict, pbc_bonds_dict = \
connect_atoms(astructure, atoms_to_connect=pair,
cov_offset=self.cov_offset, delete_existing=False,
pbc_bonding=True, only_pbc_bonds=False,
max_valencies=max_valencies)
abond = list(bonds_dict) + list(pbc_bonds_dict)
if abond:
bonds_made = True
abond = astructure.getBond(*(abond.pop()))
abond.order = pbc_bond.order
if bonds_made:
pairs = [curr_tail_head, neighbor_tail_head]
self.cleanUpPBCBonds(astructure, pairs, \
not pbc_bond.also_reg_bond)
[docs]def add_labeled_pbc_bond(astructure, atom1, atom2, order, \
is_pbc_bond=False, also_reg_bond=False, color=PBCBond.PBC_BOND_COLOR):
"""
Add the specified bond to the provided structure
and label it depending on if it is a PBC bond.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure to which the bond will be added
:type atom1: int
:param atom1: the first atom index of the bond
:type atom2: int
:param atom2: the second atom index of the bond
:type order: int
:param order: the bond order
:type is_pbc_bond: bool
:param is_pbc_bond: if True then the specified bond
is a PBC bond so we will label it as so
:type also_reg_bond: bool
:param also_reg_bond: if True indicates that the given
PBC bond is also a regular bond, meaning that one of the
atoms of the PBC bond is covalently bound to two copies of
the other atom, one inside the cell and one outside the cell
:type color: int or None
:param color: if integer specifies that PBC bonds, that are
not also regular bonds, should be colored with this color
"""
abond = msutils.add_or_update_bond_order(astructure, atom1, atom2, order)
set_representation_bond(abond)
if is_pbc_bond:
label_pbc_bond(abond, also_reg_bond=also_reg_bond, color=color)
[docs]def get_natom_btw_two_cells(cell1, cell2, extents, size):
"""
Using a traversal path of c then b then a, return the
number of atoms between two cells, of the given
size, in a super cell with the given extents.
:type cell1: tuple
:param cell1: a triple of cell indices for the first cell
:type cell2: tuple
:param cell2: a triple of cell indices for the second cell
:type extents: list
:param extents: contains the number of cells along a, b, and
c lattice vectors in the super cell
:type size: int
:param size: the number of atoms in a cell
:rtype: int
:return: the number of atoms between the two cells
"""
ncell_to_cell1 = get_collapsed_index(cell1, *extents)
ncell_to_cell2 = get_collapsed_index(cell2, *extents)
ncell_diff = ncell_to_cell2 - ncell_to_cell1
natom_diff = ncell_diff * size
return natom_diff
[docs]def delete_all_pbc_bonds(astructure):
"""
Delete all PBC bonds from the given structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure from which the bonds
will be deleted
"""
bonds = []
for abond in astructure.bond:
pbc_bond = abond.property.pop(PBCBond.PBC_BOND_KEY, None)
if pbc_bond:
also_reg_bond = abond.property.pop(PBCBond.ALSO_REG_BOND_KEY, None)
if not also_reg_bond:
bonds.append((abond.atom1.index, abond.atom2.index))
for abond in bonds:
astructure.deleteBond(*abond)
[docs]def is_pbc_bond(astructure,
atom1,
atom2,
check_also_reg_bond=False,
unit_lattice_vectors=None,
pbc=None):
"""
Return a (is_pbc_bond, also_reg_bond, bond_distance) tuple
that indicates (1) whether the specified bond is a PBC bond,
(2) if checked whether that PBC bond is also a regular bond,
and (3) the bond length.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure that has the bond
:type atom1: schrodinger.structure._StructureAtom
:param atom1: the first atom of the bond
:type atom2: schrodinger.structure._StructureAtom
:param atom2: the second atom of the bond
:type check_also_reg_bond: bool
:param check_also_reg_bond: check if the PBC bond is also
a regular bond, meaning that one of the atoms of the PBC
bond is covalently bound to two copies of the other atom,
one inside the cell and one outside the cell
:type unit_lattice_vectors: list of (numpy.array, float) tuples
or None
:param unit_lattice_vectors: contains normalized a, b, and c
lattice vectors and their original lengths, is None if no
check for PBC bonds that are also regular bonds is being
performed
:type pbc: `schrodinger.infra.structure.PBC` or None
:param pbc: The infrastructure PBC created from the Chorus box properties,
if None, pbc object will be created
:rtype: tuple
:return: a (is_pbc_bond, also_reg_bond, bond_distance) tuple
"""
threshold = PBCBond.PBC_BOND_THRESHOLD
if pbc:
pbc_distance = pbc.getDistance(astructure, atom1.index, astructure,
atom2.index)
else:
pbc_distance = pymmlibs.mmct_atom_get_distance_pbc(
astructure, atom1.index, astructure, atom2.index)
real_distance = mm.mmct_atom_get_distance(astructure.handle, atom1,
astructure.handle, atom2)
pbc_bond = also_reg_bond = False
if abs(pbc_distance - real_distance) > threshold:
pbc_bond = True
elif check_also_reg_bond:
# here we need to check if the bond in question is both
# a regular and PBC bond, for example consider the primitive
# unit cell of Si which contains two atoms that are actually
# bound to each other but which also require a PBC bond
pbc_bond_vec = numpy.array(atom2.xyz) - numpy.array(atom1.xyz)
for (avec, alen) in unit_lattice_vectors:
projection = numpy.dot(pbc_bond_vec, avec)
if abs(abs(projection) - 0.5 * alen) <= threshold:
pbc_bond = also_reg_bond = True
break
return (pbc_bond, also_reg_bond, pbc_distance)
[docs]def get_lattice_param_properties(astructure):
"""
Return a tuple containing the six lattice parameter
properties of the given structure.
:type astructure: schrodinger.Structure.structure
:param astructure: the structure that has the lattice
parameter properties defined
:rtype: tuple
:return: contains the six lattice parameter properties
"""
keys = [Crystal.A_KEY, Crystal.B_KEY, Crystal.C_KEY, \
Crystal.ALPHA_KEY, Crystal.BETA_KEY, Crystal.GAMMA_KEY]
return tuple((astructure.property[x] for x in keys))
[docs]def get_normalized_spanning_vectors(lattice_vectors):
"""
Return the thirteen unique normalized spanning vectors for
the cell.
:type lattice_vectors: list of numpy.array
:param lattice_vectors: the lattice a, b, and c vectors
:rtype: dict
:return: keys are tuples of cell index triples, values are tuples
of normalized spanning vectors and their original lengths
"""
spanning_vecs = {}
for triple in AHEAD_TRIPLES:
spanning_vec = \
sum((x[0] * x[1] for x in zip(triple, lattice_vectors)))
unit_spanning_vec = transform.get_normalized_vector(spanning_vec)
alen = transform.get_vector_magnitude(spanning_vec)
spanning_vecs[triple] = (unit_spanning_vec, alen)
return spanning_vecs
[docs]def get_element_priority(pair):
"""
Return the formula ordering priority of the element
in the given pair.
:type pair: tuple
:param pair: (element, number) tuple
:rtype: int
:return: the element priority
"""
# heavier elements are given higher priority, i.e. moved earlier
# in the formula string, elements in LAST_ELEMENTS always go last
# by heaviness
symbol = pair[0]
atomic_number = get_atomic_number(symbol)
priority = mm.MMELEMENTS_MAX - atomic_number
if symbol in LAST_ELEMENTS:
priority += mm.MMELEMENTS_MAX
return priority
[docs]def trans_cart_to_fract(cart_vec, a_param, b_param, c_param, alpha_param,
beta_param, gamma_param):
"""
Transform the given vector in the Cartesian basis to the fractional
basis.
:type cart_vec: numpy.array
:param cart_vec: the vector to be transformed from the Cartesian basis
to the fractional basis
:type a_param: float
:param a_param: the lattice a parameter
:type b_param: float
:param b_param: the lattice b parameter
:type c_param: float
:param c_param: the lattice c parameter
:type alpha_param: float
:param alpha_param: the lattice alpha parameter
:type beta_param: float
:param beta_param: the lattice beta parameter
:type gamma_param: float
:param gamma_param: the lattice gamma parameter
:rtype: numpy.array
:return: the given vector in the fractional basis
"""
# get the transform
crystal_obj = Crystal(None,
space_group=P1_SPACE_GROUP_SYMBOL,
a_param=a_param,
b_param=b_param,
c_param=c_param,
alpha_param=alpha_param,
beta_param=beta_param,
gamma_param=gamma_param)
crystal_obj.setCrystalSymmetry()
xyz2cry = numpy.array(crystal_obj.xyz2cry)
# transform the input vector
fract_vec = numpy.dot(xyz2cry, numpy.array(cart_vec))
return fract_vec
[docs]def trans_fract_to_cart(fract_vec, a_param, b_param, c_param, alpha_param,
beta_param, gamma_param):
"""
Transform the given vector in the fractional basis to the
Cartesian basis.
:type fract_vec: numpy.array
:param fract_vec: the vector to be transformed from the fractional
basis to the Cartesian basis
:type a_param: float
:param a_param: the lattice a parameter
:type b_param: float
:param b_param: the lattice b parameter
:type c_param: float
:param c_param: the lattice c parameter
:type alpha_param: float
:param alpha_param: the lattice alpha parameter
:type beta_param: float
:param beta_param: the lattice beta parameter
:type gamma_param: float
:param gamma_param: the lattice gamma parameter
:rtype: numpy.array
:return: the given vector in the Cartesian basis
"""
# get the transform
crystal_obj = Crystal(None,
space_group=P1_SPACE_GROUP_SYMBOL,
a_param=a_param,
b_param=b_param,
c_param=c_param,
alpha_param=alpha_param,
beta_param=beta_param,
gamma_param=gamma_param)
crystal_obj.setCrystalSymmetry()
cry2xyz = numpy.array(crystal_obj.cry2xyz)
# transform the input vector
cart_vec = numpy.dot(cry2xyz, numpy.array(fract_vec))
return cart_vec
[docs]def trans_cart_to_frac_from_vecs(coords, a_vec, b_vec, c_vec, rec=False):
"""
Transform coordinates from (reciprocal) Cartesian to (reciprocal) fractional
using lattice vectors.
:type coords: numpy.array
:param coords: a list of Cartesian coordinates
:type a_vec: list of 3 floats
:param a_vec: 'a' lattice vector
:type b_vec: list of 3 floats
:param b_vec: 'b' lattice vector
:type c_vec: list of 3 floats
:param c_vec: 'c' lattice vector
:type rec: bool
:param rec: If True, work in reciprocal space
:rtype: numpy array
:return: Coordinates in fractional coordinates
"""
if rec:
ra_vec, rb_vec, rc_vec = \
get_reciprocal_lattice_vectors(a_vec, b_vec, c_vec)
cry2xyz, xyz2cry = get_conv_from_vecs(ra_vec, rb_vec, rc_vec)
else:
cry2xyz, xyz2cry = get_conv_from_vecs(a_vec, b_vec, c_vec)
return numpy.dot(coords, xyz2cry)
[docs]def trans_frac_to_cart_from_vecs(coords, a_vec, b_vec, c_vec, rec=False):
"""
Transform coordinates from (reciprocal) fractional to (reciprocal) Cartesian
using lattice vectors.
:type coords: numpy.array
:param coords: a list of fractional coordinates
:type a_vec: list of 3 floats
:param a_vec: 'a' lattice vector
:type b_vec: list of 3 floats
:param b_vec: 'b' lattice vector
:type c_vec: list of 3 floats
:param c_vec: 'c' lattice vector
:type rec: bool
:param rec: If True, work in reciprocal space
:rtype: numpy array
:return: Coordinates in Cartesian coordinates
"""
if rec:
ra_vec, rb_vec, rc_vec = \
get_reciprocal_lattice_vectors(a_vec, b_vec, c_vec)
cry2xyz, xyz2cry = get_conv_from_vecs(ra_vec, rb_vec, rc_vec)
else:
cry2xyz, xyz2cry = get_conv_from_vecs(a_vec, b_vec, c_vec)
return numpy.dot(coords, cry2xyz)
[docs]def get_cell(asu,
space_group=None,
lattice_params=None,
extents=None,
xtal_kwargs=None):
"""
Build and return a crystalline cell using the given
asymmetric unit, space group, lattice parameters, and
extents.
:type asu: `schrodinger.structure.Structure`
:param asu: the asymmetric unit
:type space_group: str or None
:param space_group: the full or short Hermann-Mauguin symbol
of the space group or None in which case the asu structure
property Crystal.SPACE_GROUP_KEY will be used
:type lattice_params: list or None
:param lattice_params: the six lattice parameters or None in
which case the asu structure properties Crystal.A_KEY,
Crystal.B_KEY, Crystal.C_KEY, Crystal.ALPHA_KEY, Crystal.BETA_KEY,
and Crystal.GAMMA_KEY will be used
:type extents: list or None
:param extents: the integer extents along the a, b, and c lattice
vectors or None if there are none
:type xtal_kwargs: dict or None
:param xtal_kwargs: extra xtal.Crystal kwargs or None if there
are none
:rtype: `schrodinger.structure.Structure`
:return: the built cell
"""
# initialize kwargs
if xtal_kwargs is not None:
kwargs = dict(xtal_kwargs)
else:
kwargs = {}
# add space group
if space_group is not None:
kwargs['space_group'] = space_group
# add lattice parameters
if lattice_params is not None:
kwargs.update(dict(zip(LATTICE_PARAMS_KEYS, lattice_params)))
# add extents
if extents is not None:
kwargs.update(dict(zip(EXTENTS_KEYS, extents)))
acrystal = Crystal(asu, **kwargs)
acrystal.orchestrate()
return acrystal.crystal_super_cell
[docs]def get_vectors_from_chorus(astructure):
"""
Return the three lattice vectors from the nine
lattice chorus properties.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure that has the chorus
properties
:raise ValueError: if any chorus property is missing
:rtype: numpy.array
:return: the three lattice vectors
"""
afunc = lambda x: numpy.array([astructure.property[key] for key in x])
try:
vectors = list(map(afunc, Crystal.CHORUS_BOX_KEY_VECTORS))
except KeyError:
err = ('The given structure is missing one or more '
'lattice chorus properties.')
raise ValueError(err)
return numpy.array(vectors)
[docs]def get_conv_from_vecs(a_vec, b_vec, c_vec):
"""
Generate matrices to convert from fractional to Cartesian and back.
:type a_vec: list of 3 floats
:param a_vec: 'a' lattice vector
:type b_vec: list of 3 floats
:param b_vec: 'b' lattice vector
:type c_vec: list of 3 floats
:param c_vec: 'c' lattice vector
:rtype: tuple of matrices
:return: Matrices that convert fractional atomic coordinates to Cartesian
and back
"""
cry2xyz = numpy.empty(shape=[3, 3])
cry2xyz[0] = a_vec.copy()
cry2xyz[1] = b_vec.copy()
cry2xyz[2] = c_vec.copy()
xyz2cry = numpy.linalg.inv(cry2xyz)
return cry2xyz, xyz2cry
[docs]def create_new_box(struct, minval=0., buffer=PBC_BUFFER_LENGTH):
"""
Create a new box that is large enough to encompass the X, Y and Z
lengths of the system
:type struct: `schrodinger.structure.Structure`
:param struct: The structure to work on
:type minval: float
:param minval: Minimum length of the cell edge or 0. if vdW is fine
:param float buffer: The buffer to add to all PBC lengths. We add 2.0 as the
default vdw buffer to make sure atoms at the boundary don't clash with
their mirror images.
"""
assert buffer >= PBC_BUFFER_LENGTH, (
"A minimum of 2.0 (A) vdw buffer is required so atoms "
"at the boundary don't clash with their mirror images")
vals = struct.getXYZ().ptp(0)
vals = [max(val + buffer, minval) for val in vals]
store_chorus_box_props(struct, vals[0], by=vals[1], cz=vals[2])
[docs]def set_lattice_properties(astructure, lattice_properties):
"""
Set the given lattice properties on the structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure on which to set the lattice
properties
:type lattice_properties: list
:param lattice_properties: a, b, c, alpha, beta, and gamma
lattice properties
"""
keys = [Crystal.A_KEY, Crystal.B_KEY, Crystal.C_KEY, \
Crystal.ALPHA_KEY, Crystal.BETA_KEY, Crystal.GAMMA_KEY]
for key, prop in zip(keys, lattice_properties):
astructure.property[key] = prop
[docs]def make_p1(astructure, logger=None, in_place=False):
"""
Make a P1 cell.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure to make P1
:type logger: `logging.Logger` or None
:param logger: if not None then the logger for printing
:type in_place: bool
:param in_place: if True then operate directly on the given
structure as opposed to a copy of it
:rtype: `schrodinger.structure.Structure`
:return: the P1 cell
"""
if in_place:
in_astructure = astructure
else:
in_astructure = astructure.copy()
try:
chorus_properties = get_chorus_properties(in_astructure)
except KeyError:
chorus_properties = None
if chorus_properties is None:
try:
lattice_properties = get_lattice_param_properties(in_astructure)
except KeyError:
create_new_box(in_astructure)
chorus_properties = get_chorus_properties(in_astructure)
else:
chorus_properties = get_chorus_from_params(lattice_properties)
lattice_properties = get_params_from_chorus(chorus_properties)
set_lattice_properties(in_astructure, lattice_properties)
clear_asu_and_fractional_properties(in_astructure)
in_astructure.property[Crystal.SPACE_GROUP_KEY] = P1_SPACE_GROUP_SYMBOL
in_astructure.property[Crystal.SPACE_GROUP_ID_KEY] = P1_SPACE_GROUP_ID
volume_key = Crystal.UNIT_CELL_VOLUME_KEY
if volume_key in in_astructure.property:
volume = get_volume_from_params(lattice_properties)
in_astructure.property[volume_key] = volume
formula_key = Crystal.UNIT_CELL_FORMULA_KEY
if formula_key in in_astructure.property:
formula = get_unit_cell_formula(in_astructure)
in_astructure.property[formula_key] = formula
return in_astructure
[docs]def sync_pbc(st, create_pbc=False, in_place=False):
"""
Return the given structure with a synchronized PBC, if create_pbc
is True and the structure lacks a PBC then one will be created,
otherwise this function will return False if there is no PBC.
:type st: structure.Structure
:param st: the structure with the PBC to be synchronized
:type create_pbc: bool
:param create_pbc: if True and the given structure lacks a PBC
then a minimal PBC will be created, if False and if the given
structure lacks a PBC then this function will return False
:type in_place: bool
:param in_place: if True then operate directly on the given
structure as opposed to a copy of it
:rtype: structure.Structure or bool
:return: the structure with a synchronized PBC or False if there
was no PBC and one wasn't created
"""
try:
lattice_params = get_lattice_param_properties(st)
except KeyError:
lattice_params = None
try:
chorus_params = get_chorus_properties(st)
except KeyError:
chorus_params = None
if lattice_params or chorus_params or create_pbc:
st = make_p1(st, in_place=in_place)
lattice_params = get_lattice_param_properties(st)
chorus_params = get_chorus_from_params(lattice_params)
set_pbc_properties(st, chorus_params)
return st
else:
return False
[docs]def sync_pbc2(struct,
lattice_params=None,
chorus_params=None,
prioritize_cparams=True):
"""
Sync PBC properties in place (without creating a new structure) for struct.
If all PBC properties are absent (both chorus and PDB) return False. It is
possible to provide new lattice or chorus parameters and to prioritize one
of the sets.
:type: `schrodinger.structure.Structure`
:param: Structure to be modified
:type lattice_params: list or numpy.array or None
:param lattice_params: contains the a, b, c, alpha, beta, and gamma lattice
parameters or None. These will be used instead of ones possibly obtained
from the struct
:type chorus_params: list or numpy.array or None
:param chorus_params: contains the nine chorus properties,
i.e. ax, ay, az, bx, ..., cz or None. These will be used instead of ones
possibly obtained from the struct
:type prioritize_cparams: bool
:param: Prioritize chorus params over lattice params if True. If False,
lattice params have priority
:rtype: bool
:return: True on syncing success, False if both sets were not provided or
empty
"""
try:
if lattice_params is None:
lattice_params = get_lattice_param_properties(struct)
except KeyError:
lattice_params = []
try:
if chorus_params is None:
chorus_params = get_chorus_properties(struct)
except KeyError:
chorus_params = []
if not any(map(len, [lattice_params, chorus_params])):
return False
use_cparams = prioritize_cparams and len(chorus_params)
if not use_cparams and len(lattice_params):
pbc = infrastructure.PBC(*lattice_params)
vecs = numpy.array([v for v in pbc.getBoxVectors()])
set_pbc_properties(struct, vecs.flat)
else:
set_pbc_properties(struct, chorus_params)
return True
[docs]def get_pbc_bonds_dict(astructure):
"""
Return a PBC bonds dict for the given structure.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure with the PBC bonds from
which to create the dict
:rtype: dict
:return: keys are tuples of PBC bond atom index pairs,
values are tuples of bond orders and booleans indicating whether
the PBC bond is also a regular bond
"""
pbc_bonds = {}
for abond in astructure.bond:
if abond.property.get(PBCBond.PBC_BOND_KEY, None):
pair = (abond.atom1.index, abond.atom2.index)
also_reg = abond.property.get(PBCBond.ALSO_REG_BOND_KEY, None)
pbc_bonds[pair] = (abond.order, also_reg)
return pbc_bonds
[docs]def set_representation_bond(abond):
"""
Set the representation of the given bond.
:type abond: schrodinger.structure._StructureBond
:param abond: the bond to set the representation for
"""
# see MAE-36614 where it was decided to use representation of the
# first atom in the bond to set the representation of the bond
abond.setStyle(ATOM_TO_BOND_REPR_DICT[abond.atom1.style])
[docs]def set_pbc_properties(astructure, chorus_properties):
"""
Set the chorus and PDB properties on the given structure using
the given chorus properties.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure on which the properties are set
:type chorus_properties: list
:param chorus_properties: contains the nine chorus properties,
i.e. ax, ay, az, bx, ..., cz
"""
if len(chorus_properties) != 9:
raise ValueError(
"Expected 9 chorus box values, found {}".format(chorus_properties))
pbc = infrastructure.PBC(*chorus_properties)
pbc.applyToStructure(astructure)
[docs]def get_simple_supercell_matrix(supercell_matrix):
"""
Get minimal diagonal elements of a simple supercell matrix starting from
(non)-diagonal supercell matrix. Supercell matrix can be:
-1, 1, 1
2, -3, 4
-5, 6, 7
Resulting simple supercell (described by the matrix) must be able to hold
the supercell above.
In the case of a smaller lattice, for example:
0.5 0 0
0 0.5 0
0 0 0.5
[1, 1, 1] will be returned
:type supercell_matrix: 3 x 3 list of floats
:param supercell_matrix: Supercell matrix
:rtype: list of three ints
:return: Diagonal elements of the simple supercell matrix
(must be positive nonzero integers)
"""
matrix = numpy.array(supercell_matrix)
if numpy.linalg.det(matrix) <= 0.0:
raise ValueError('Determinant of the supercell matrix must have a '
'nonzero positive value.')
verts = numpy.array([[0.0, 0.0, 0.0]])
for indx in range(3):
verts = numpy.vstack([verts, matrix[indx]])
# Vector sum of unit cell vectors is computed to get "far" vertices
verts = numpy.vstack([verts, matrix[0] + matrix[1]])
verts = numpy.vstack([verts, matrix[1] + matrix[2]])
verts = numpy.vstack([verts, matrix[0] + matrix[2]])
verts = numpy.vstack([verts, matrix[0] + matrix[1] + matrix[2]])
ret = verts.ptp(axis=0).astype(int)
ret = [x + 1 if x < 1 else x for x in ret]
return ret
[docs]def get_simple_supercell(struct, supercell_matrix):
"""
Get supercell structure.
:type struct: `structure.Structure`
:param struct: Input structure
:type supercell_matrix: List of 3 floats
:param supercell_matrix: Diagonal elements of the supercell matrix
:rtype: `structure.Structure`
:return: Supercell structure
"""
ret_struct = struct.copy()
xyz = struct.getXYZ()
vecs = numpy.array(get_vectors_from_chorus(struct))
new_vecs = numpy.dot(numpy.diag(supercell_matrix), vecs)
set_pbc_properties(ret_struct, new_vecs.flat)
for indexa in range(supercell_matrix[0]):
for indexb in range(supercell_matrix[1]):
for indexc in range(supercell_matrix[2]):
if not any([indexa, indexb, indexc]):
continue
new_vecs = numpy.dot(numpy.diag([indexa, indexb, indexc]), vecs)
translate_vec = new_vecs.sum(axis=0).tolist()
next_cell = struct.copy()
next_cell.setXYZ(xyz + translate_vec)
ret_struct.extend(next_cell)
return ret_struct
[docs]def get_physical_properties(struct, vecs=None):
"""
Get cell formula, volume and density of a struct.
:type: `structure.Structure`
:param: Input structure
:type vecs: list(list) or None
:param vecs: Lattice vectors. If None, will be obtained from structure
:rtype: float, float, float
:return: formula, volume and density
"""
formula = get_unit_cell_formula(struct)
if vecs is None:
vecs = get_vectors_from_chorus(struct)
volume = get_volume_from_vecs(vecs)
# do not use the structure.total_weight or atom.atomic_weight attrs as they
# may be wrong given that our structure may have open valencies and those
# attrs include implicit hydrogens, in fact atom.atomic_weight is safe for
# coarse grain particles
volume_cm3 = volume * math.pow(1.0 / Crystal.CM_TO_ANGSTROM, 3)
if not msutils.is_coarse_grain(struct):
mass = sum(
[get_atomic_weight(atom.atomic_number) for atom in struct.atom])
else:
mass = sum([atom.atomic_weight for atom in struct.atom])
density = mass / constants.N_A / volume_cm3
return formula, volume, density
[docs]def set_physical_properties(struct):
"""
Set cell formula, volume and density to struct.
:type: `structure.Structure`
:param: Input structure
"""
formula, volume, density = get_physical_properties(struct)
struct.property[Crystal.UNIT_CELL_FORMULA_KEY] = formula
struct.property[Crystal.UNIT_CELL_VOLUME_KEY] = volume
struct.property[Crystal.UNIT_CELL_DENSITY_KEY] = density
[docs]def pdist_vec_row_col(d, i):
"""
Convert from triangular indices of distance matrix to indices of the square
form.
:type d: int
:param d: row length of the original triangular matrix
:type i: numpy.array(int)
:param i: Condensed triangular indices, 0-indexed
:rtype: numpy.array(int), numpy.array(int)
:return: row and column indices (0-indexed) from the corresponding square
form
"""
# Based on https://stackoverflow.com/q/5323818/
b = 1. - 2. * d
x = numpy.floor((-b - numpy.sqrt(b**2. - 8. * i)) / 2.).astype(int)
y = (i + x * (b + x + 2.) / 2. + 1.).astype(int)
return x, y
[docs]def preserve_bonds(struct, to_keep, to_remove):
"""
Try to preserve bonding during the remove of the duplicates atoms.
:type struct: structure.Structure
:param struct: Structure to preserve bonding for
:type to_keep: list
:param to_keep: List of atom indices to keep
:type to_remove: list(list)
:param to_remove: List of lists of atom indices to be removed. Each item in
the list is a list of duplicated atom indices corresponding to the atom
index from to_keep list
"""
# This function gives wrong results for the case (see RB of MATSCI-6684):
# When the only bond instance that exists in the structures
# involves pairs of duplicate atoms being removed, for example a bond
# between a duplicate of atom 1 and a duplicate of atom 2 with no other
# instance of the bond in existence.
assert len(to_keep) == len(to_remove)
if len(to_keep) == 0:
return
to_remove_set = set(numpy.concatenate(to_remove))
# Delete bonds between atoms from to_remove and all the rest to-be-deleted
# atoms
bonds_to_delete = []
for remove_jdxs in to_remove:
for jdx in remove_jdxs:
atom2 = struct.atom[jdx]
bonds_to_delete.extend([
frozenset((atom2, batom))
for batom in atom2.bonded_atoms
if batom.index in to_remove_set
])
build.delete_bonds(list(set(bonds_to_delete)))
max_valencies = elementalprops.get_mmct_max_valencies()
for keep_idx, remove_jdxs in zip(to_keep, to_remove):
atom1 = struct.atom[keep_idx]
# Delete bonds between atom1 and all to-be-deleted atoms
build.delete_bonds([(atom1, batom)
for batom in atom1.bonded_atoms
if batom.index in to_remove_set])
for jdx in remove_jdxs:
atom2 = struct.atom[jdx]
# to_remove atom2 is only bonded to the to_keep atom1 at this point
# bonds may be removed inside the loop, make list of bonded atoms
batoms = list(atom2.bonded_atoms)
for batom in batoms:
# batom is one of the atoms to be kept
bond = struct.getBond(atom2, batom)
# Save bond type and styles before deleting it
bond_type = bond.type
to_style, from_style = bond.to_style, bond.from_style
struct.deleteBond(atom2, batom)
# Prevent overbonding
if (atom1.bond_total < max_valencies[atom1.element] and
batom.bond_total < max_valencies[batom.element]):
bond = msutils.add_or_update_bond_type(
struct, atom1, batom, bond_type)
bond.to_style, bond.from_style = to_style, from_style
[docs]def translate_to_1st_cell(coords,
vecs,
is_cartesian=True,
fract_offset=FRACT_OFFSET,
axes=None):
"""
Translate coordinates to the 1st cell (in range [0, 1). If it is a
supercell with original vectors, overlap atoms are expected.
:param list coords: Coordinates (can be Cartesian or fractional,
see is_cartesian)
:param list vecs: Lattice vectors
:param bool is_cartesian: Whether input coordinates are Cartesian or
fractional
:param float fract_offset: The threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:type axes: iterable or None
:param axes: List of axes to use
:return: Translated Cartesian coordinates
"""
# Currently it translates in range [0 - fract_offset, 1) (MATSCI-6647)
if is_cartesian:
coords = trans_cart_to_frac_from_vecs(coords, *vecs)
# Move into range (-1, 1)
new_fracs = numpy.fmod(coords, 1.)
# Move into range [0, 1]
new_fracs[new_fracs < 0.] += 1.
# Enforce range [0, 1)
new_fracs[new_fracs >= 1 - fract_offset] -= 1.
tmp_fracs = numpy.array(coords)
tmp_fracs[:, axes] = new_fracs[:, axes]
new_fracs = tmp_fracs
return trans_frac_to_cart_from_vecs(new_fracs, *vecs)
[docs]def coords_outside_1st_cell(coords,
vecs,
is_cartesian=True,
fract_offset=FRACT_OFFSET,
axes=None):
"""
Returns coordinate indices (atom indices - 1) of the atoms outside of the
cell defined with vecs.
:param list coords: Coordinates (can be Cartesian or fractional,
see is_cartesian)
:param list vecs: Lattice vectors
:param bool is_cartesian: Whether input coordinates are Cartesian or
fractional
:param float fract_offset: The threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:type axes: iterable or None
:param axes: List of axes to use
:rtype: numpy.array
:return: Indices outside the first cell
"""
if is_cartesian:
coords = trans_cart_to_frac_from_vecs(coords, *vecs)
# Atom indices outside the first cell
indices = numpy.where(
numpy.logical_or(coords[:, axes] < 0. - fract_offset,
coords[:, axes] >= 1.))[0]
return indices
[docs]def is_normal_surface(vecs):
"""
Checks if C-axis is normal to the a-b plane from lattice vectors.
:param numpy.array vecs: Lattice vectors
:rtype: bool
:return: Whether C-axis is normal to the a-b plane
"""
params = get_params_from_chorus(vecs.flat)
return numpy.isclose(params[3:5], [90., 90.]).all()
[docs]def move_atoms_into_cell(struct,
frac_coords=None,
overlap_tresh=OVERLAP_ATOM_THRESHOLD,
fract_offset=FRACT_OFFSET,
preserve_bonding=False):
"""
Get structure with all the atoms moved into the first cell.
:type: `structure.Structure`
:param: Input structure
:type frac_coords: numpy arrays of arrays of 3 floats or None
:param frac_coords: Fractional coordinates
:type overlap_tresh: float
:param overlap_tresh: Distance between two atoms, such that they are
considered overlapping
:type fract_offset: float
:param fract_offset: The threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
:type preserve_bonding: bool
:param preserve_bonding: If True, preserve bonding between atoms
(might be slow)
:rtype: `structure.Structure`
:return: Structure with all the atoms moved inside first unit cell
"""
vecs = get_vectors_from_chorus(struct)
if frac_coords is not None:
if len(frac_coords) != struct.atom_total:
raise ValueError('Number of fractional coordinates (%d) is '
'different than number of atoms (%d).' %
(len(frac_coords), struct.atom_total))
else:
frac_coords = trans_cart_to_frac_from_vecs(struct.getXYZ(), *vecs)
# Copy original structure
new_struct = struct.copy()
new_carts = translate_to_1st_cell(frac_coords,
vecs,
is_cartesian=False,
fract_offset=fract_offset)
new_struct.setXYZ(new_carts)
# Before scipy.distance.pdist was used which places entire distance matrix
# in the memory at once (MATSCI-7197)
pbc = infrastructure.PBC(*get_chorus_properties(struct))
atoms_to_keep, atoms_to_delete = get_duplicate_atoms(
new_struct, pbc=pbc, duplicate_thresh=overlap_tresh)
if len(atoms_to_delete):
if preserve_bonding:
preserve_bonds(new_struct, atoms_to_keep, atoms_to_delete)
new_struct.deleteAtoms(set(numpy.concatenate(atoms_to_delete)))
label_pbc_bonds(new_struct)
clear_asu_and_fractional_properties(new_struct)
set_physical_properties(new_struct)
return new_struct
[docs]def get_unit_lattice_vector_info(astructure):
"""
Return a list of tuples containing unit lattice vector information,
i.e. the unit lattice vectors and their original lengths.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure with lattice vectors defined by
chorus box properties
:rtype: list of tuples
:return: each (numpy.array, float) tuple contains (1) the unit
lattice vector and (2) the length of the original vector
"""
lattice_vectors = [
numpy.array([astructure.property[key]
for key in x])
for x in Crystal.CHORUS_BOX_KEY_VECTORS
]
unit_lattice_vectors = \
[(transform.get_normalized_vector(avec), \
transform.get_vector_magnitude(avec)) for avec in lattice_vectors]
return unit_lattice_vectors
[docs]def label_pbc_bonds(astructure, pbc=None, check_also_reg_bond=None, is_cg=None):
"""
Label PBC bonds.
:type astructure: `schrodinger.structure.Structure`
:param astructure: the structure with the bonds to label
:type pbc: `schrodinger.infra.structure.PBC` or None
:param pbc: The infrastructure PBC created from the Chorus box properties,
if None, pbc object will be created
:type check_also_reg_bond: bool or None
:param check_also_reg_bond: if True then PBC bonds will be
checked to see if they are also regular bonds, meaning that
one of the atoms of the PBC bond is covalently bound to two
copies of the other atom, one inside the cell and one outside
the cell. If None, this will be evaluated
:type is_cg: bool or None
:param is_cg: True, if structure is CG, otherwise False. If None, structure
will be evaluated (slow)
:type pbc: `schrodinger.infra.structure.PBC` or None
:param pbc: The infrastructure PBC created from the Chorus box properties,
if None, pbc object will be created
"""
unit_lattice_vectors = get_unit_lattice_vector_info(astructure)
if not pbc:
pbc = infrastructure.PBC(*get_chorus_properties(astructure))
if check_also_reg_bond is None:
if is_cg is None:
is_cg = msutils.is_coarse_grain(astructure)
check_also_reg_bond = get_check_also_reg_bond(astructure,
pbc=pbc,
is_cg=is_cg)
for abond in astructure.bond:
abond.property.pop(PBCBond.PBC_BOND_KEY, None)
abond.property.pop(PBCBond.ALSO_REG_BOND_KEY, None)
pbc_bond, also_reg_bond, distance = is_pbc_bond(
astructure,
abond.atom1,
abond.atom2,
check_also_reg_bond=check_also_reg_bond,
unit_lattice_vectors=unit_lattice_vectors,
pbc=pbc)
if pbc_bond:
label_pbc_bond(abond, also_reg_bond=also_reg_bond, color=None)
[docs]def label_pbc_bond(abond, also_reg_bond=False, color=PBCBond.PBC_BOND_COLOR):
"""
Label this PBC bond.
:type abond: schrodinger.structure._StructureBond
:param abond: the PBC bond to label
:type also_reg_bond: bool
:param also_reg_bond: if True indicates that the given
PBC bond is also a regular bond, meaning that one of the
atoms of the PBC bond is covalently bound to two copies of
the other atom, one inside the cell and one outside the cell
:type color: int or None
:param color: if integer specifies that a PBC bond, that is
not also a regular bond, should be colored with this color,
if None then no coloring is performed
"""
abond.property[PBCBond.PBC_BOND_KEY] = True
if also_reg_bond:
abond.property[PBCBond.ALSO_REG_BOND_KEY] = True
elif color is not None:
abond.property[PBCBond.PBC_BOND_COLOR_KEY] = color
[docs]def get_primitive_cell(struct_in, set_space_group=False):
"""
Get primitive cell.
:type struct: `schrodinger.structure.Structure`
:param struct: Original structure
:type set_space_group: bool
:type set_space_group: If True, set space group in the structure property
:rtype: `schrodinger.structure.Structure`
:return: Primitive cell with the space group set, if requested
"""
vecs = numpy.array(get_vectors_from_chorus(struct_in))
fcoords = trans_cart_to_frac_from_vecs(struct_in.getXYZ(), *vecs)
anums = [a.atomic_number for a in struct_in.atom]
spg_cell = (vecs, fcoords, anums)
prim_cell = spglib.find_primitive(spg_cell, symprec=ASSIGN_SPG_SYMPREC)
struct = struct_in.copy()
if prim_cell is None:
return struct
struct.deleteAtoms(range(1, struct.atom_total + 1))
vecs, fcoords, anums = prim_cell
lparams = get_params_from_chorus(vecs.flat)
vecs = numpy.array(get_lattice_vectors(*lparams))
set_pbc_properties(struct, vecs.flat)
coords = trans_frac_to_cart_from_vecs(fcoords, *vecs)
for xyz, anum in zip(coords, anums):
struct.addAtom(msutils.get_atomic_element(anum), *xyz)
set_physical_properties(struct)
struct.property[PBC_POSITION_KEY] = ANCHOR_PBC_POSITION % ('0', '0', '0')
if set_space_group:
assign_space_group(struct)
return struct
[docs]def get_std_cell_from_spglib_dataset(struct, dataset):
"""
Get standardized cell (structure) from spglib dataset and copy all the
structure properties (except symmetry related) from struct_in.
:type struct: `schrodinger.structure.Structure`
:param struct: Structure to get structure properties from
:type dataset: dict
:param dataset: Dataset containing cell and space group properties, see:
https://atztogo.github.io/spglib/python-spglib.html#get-symmetry-dataset
:rtype: `schrodinger.structure.Structure`
:return: Standardized cell with the space group set
"""
std_struct = struct.copy()
std_struct.deleteAtoms(range(1, struct.atom_total + 1))
std_fcoords = dataset['std_positions']
std_vecs = dataset['std_lattice']
set_pbc_properties(std_struct, std_vecs.flat)
std_coords = trans_frac_to_cart_from_vecs(std_fcoords, *std_vecs)
for xyz, anum in zip(std_coords, dataset['std_types']):
std_struct.addAtom(msutils.get_atomic_element(anum), *xyz)
set_physical_properties(std_struct)
std_struct.property[PBC_POSITION_KEY] = \
ANCHOR_PBC_POSITION % ('0', '0', '0')
assign_space_group(std_struct)
return std_struct
[docs]def assign_space_group(struct_in,
symprec=ASSIGN_SPG_SYMPREC,
search_alt_cells=False):
"""
Set space group and space group id to the input structure. An error message
is returned on failure.
:type: `schrodinger.structure.Structure`
:param: Input structure that will be modified
:type symprec: float
:param symprec: Symmetry tolerance used for atomic coordinates to assign
space group
:type search_alt_cells: bool
:param search_alt_cells: If True, search for alternative cells
:rtype: list
:return: If search_alt_cells is True, conventional and primitive cells are
returned. Otherwise, empty list.
"""
def is_supercell(dataset):
eps = 1e-5 # spglib spits not precise numbers
return numpy.linalg.det(dataset['transformation_matrix']) > 1. + eps
def get_spglib_cell(struct):
struct = move_atoms_into_cell(struct)
vecs = get_vectors_from_chorus(struct)
fracs = trans_cart_to_frac_from_vecs(struct.getXYZ(), *vecs)
anums = [a.atomic_number for a in struct.atom]
return vecs, fracs, anums
spg_cell = get_spglib_cell(struct_in)
dataset = spglib.get_symmetry_dataset(spg_cell, symprec=symprec)
if not dataset:
struct_in.property[SPACE_GROUP_KEY] = P1_SPACE_GROUP_SYMBOL
struct_in.property[SPACE_GROUP_ID_KEY] = P1_SPACE_GROUP_ID
return [struct_in.copy(), struct_in.copy()] if search_alt_cells else []
if is_supercell(dataset):
# Deal with super cell. Get conventional cell and use its space group
struct = get_std_cell_from_spglib_dataset(struct_in, dataset)
spg_cell = get_spglib_cell(struct)
dataset = spglib.get_symmetry_dataset(spg_cell, symprec=symprec)
spg_id = int(dataset['number'])
symmops = space_groups.get_symmops_from_spglib(dataset)
spg_obj = None
for spg in space_groups.get_spacegroups().spg_objs[spg_id]:
if space_groups.equal_rotations(spg.symmetry_opers, symmops):
spg_obj = spg
break
if spg_obj:
struct_in.property[SPACE_GROUP_KEY] = spg_obj.space_group_short_name
struct_in.property[SPACE_GROUP_ID_KEY] = spg_obj.space_group_id
else:
struct_in.property[SPACE_GROUP_KEY] = P1_SPACE_GROUP_SYMBOL
struct_in.property[SPACE_GROUP_ID_KEY] = P1_SPACE_GROUP_ID
if not search_alt_cells:
return []
conv_cell = get_std_cell_from_spglib_dataset(struct_in, dataset)
conv_cell.title += '_conventional'
prim_cell = get_primitive_cell(struct_in, set_space_group=True)
prim_cell.title += '_primitive'
return conv_cell, prim_cell
[docs]def get_normal_surf(struct,
restore_vacuum=True,
overlap_threshold=OVERLAP_ATOM_THRESHOLD):
"""
Enforce system to have C-axis normal to the a-b plane. This is done by
straining the structure.
:type struct: `schrodinger.structure.Structure`
:param struct: Input structure
:type restore_vacuum: bool
:param restore_vacuum: If True, set vacuum value in Z direction to the original
value, before the surface is "normalized". Should be set to True for
infinite systems. If False, don't modify cell after normalization
:param float overlap_threshold: (Ang) distance used to define overlapping
atoms
:rtype: `schrodinger.structure.Structure`
:return: Output structure with C-axis normal to the a-b plane
"""
vecs = numpy.array(get_vectors_from_chorus(struct))
params = get_params_from_chorus(vecs.flat)
if numpy.isclose(params[3:5], [90., 90.]).all():
# c-axis is already normal to a-b plane
return struct.copy()
# Get the original amount of vacuum
xmax, ymax, zmax = struct.getXYZ().max(axis=0)
c_edge_coords = trans_frac_to_cart_from_vecs(numpy.array([
0.,
0.,
1.,
]), *vecs)
vacuum = c_edge_coords[2] - zmax
# This (alpha = beta = 90 ) ensures that a-b plane is normal to
# the c-axis
params[3:5] = [90., 90.]
ortho_vecs = numpy.array(get_chorus_from_params(params)).reshape(3, 3)
ortho_tmatrix = get_transformation_matrix(vecs, ortho_vecs)
# cell_only mode was proposed in MATSCI-5181
new_struct = transform_pbc(struct,
ortho_tmatrix,
cell_only=True,
overlap_threshold=overlap_threshold)
# Return here, since below vacuum is changed
if not restore_vacuum:
return new_struct
# Remove vacuum, set c-axis length to max Z coordinate of the
# atom
xmax, ymax, zmax = new_struct.getXYZ().max(axis=0)
ortho_params = list(get_lattice_param_properties(new_struct))
# Set vacuum to the original amount to have correct c-axis extension
ortho_params[2] = zmax + vacuum
# Sync chorus with PDB lattice params
chorus_params = get_chorus_from_params(ortho_params)
set_pbc_properties(new_struct, chorus_params)
clear_asu_and_fractional_properties(new_struct)
set_physical_properties(new_struct)
return new_struct
[docs]def get_normal_surf_from_HKL(hkl, vecs, max_normal_search=10):
"""
Calculate transformation matrix with c-axis most normal to a-b plane from
HKL plane indices and lattice vectors.
:type hkl: list of 3 integers
:param hkl: Miller plane indices
:type vecs: list of 3 lists of floats
:param vecs: Lattice vectors
:type max_normal_search: int
:param max_normal_search: Maximum number of linear combinations of lattice
vectors to be considered when search most normal surface
:rtype: bool, numpy.array 3 x 3 of integers
:return: True, if normal is found, otherwise False and transformation matrix
"""
# This is based on pymatgen/core/surface.py :: SlabGenerator (MIT license)
# Calculate the surface normal using the reciprocal lattice vector.
vec_hkl = trans_frac_to_cart_from_vecs(hkl, *vecs, rec=True)
normal = transform.get_normalized_vector(vec_hkl)
found_normal = False
slab_scale_factor = []
non_orth_ind = []
eye = numpy.eye(3, dtype=numpy.int)
for idx, jdx in enumerate(hkl):
if jdx == 0:
# Lattice vector is perpendicular to surface normal, i.e.,
# in plane of surface. We will simply choose this lattice
# vector as one of the basis vectors.
slab_scale_factor.append(eye[idx])
else:
#Calculate projection of lattice vector onto surface normal.
vec_length = transform.get_vector_magnitude(vecs[idx])
d_length = old_div(abs(numpy.dot(normal, vecs[idx])), vec_length)
non_orth_ind.append((idx, d_length))
# We want the vector that has maximum magnitude in the
# direction of the surface normal as the c-direction.
# Results in a more "orthogonal" unit cell.
c_index, dist = max(non_orth_ind, key=lambda t: t[1])
if len(non_orth_ind) > 1:
lcm_miller = lcm([hkl[i] for i, d in non_orth_ind])
for (idx, didx), (jdx, djdx) in itertools.combinations(non_orth_ind, 2):
vec = [0, 0, 0]
vec[idx] = -int(round(old_div(lcm_miller, hkl[idx])))
vec[jdx] = int(round(old_div(lcm_miller, hkl[jdx])))
slab_scale_factor.append(vec)
if len(slab_scale_factor) == 2:
break
index_range = sorted(reversed(
list(range(-max_normal_search, max_normal_search + 1))),
key=lambda x: abs(x))
candidates = []
for uvw in itertools.product(index_range, index_range, index_range):
if not any(uvw):
continue
# see MATSCI-4319 - where numpy.linalg.det raises RuntimeWarning
# any time the determinant is zero which is properly handled by
# this code, related to https://github.com/numpy/numpy/issues/8529
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
normal_surf_hkl = abs(numpy.linalg.det(slab_scale_factor + [uvw]))
if normal_surf_hkl < NORMAL_SURF_HKL_THR:
continue
vec = trans_cart_to_frac_from_vecs(uvw, *vecs)
vec_length = transform.get_vector_magnitude(vec)
cosine = abs(old_div(numpy.dot(vec, normal), vec_length))
candidates.append((uvw, cosine, vec_length))
if abs(abs(cosine) - 1) < NORMAL_SURF_HKL_THR:
# If cosine of 1 is found, no need to search further.
found_normal = True
break
# We want the indices with the maximum absolute cosine,
# but smallest possible length.
uvw, cosine, vec_length = max(candidates, key=lambda x: (x[1], -x[2]))
slab_scale_factor.append(uvw)
slab_scale_factor = numpy.array(slab_scale_factor)
# Let's make sure we have a left-handed crystallographic system
if numpy.linalg.det(slab_scale_factor) < 0:
slab_scale_factor *= -1
# Make sure the slab_scale_factor is reduced to avoid
# unnecessarily large slabs
reduced_scale_factor = [reduce_vector(v) for v in slab_scale_factor]
slab_scale_factor = numpy.array(reduced_scale_factor)
return found_normal, numpy.array(slab_scale_factor).T
[docs]def lcm(numbers):
"""
Get the lowest common multiple of a sequence of numbers.
:type numbers: list of integers
:param numbers: Sequence of integers
:rtype: int
:return: Lowest common multiple
"""
ret = 1
for idx in numbers:
ret = (idx * ret) // gcd(idx, ret)
return ret
[docs]def reduce_vector(vector):
"""
Get reduced vector of a transformation matrix.
:type vector: list of integers
:param vector: Components of a vector of transformation matrix
:rtype: list of integers
:return: Components of a reduced vector of transformation matrix
"""
d = abs(reduce(gcd, vector))
vector = tuple([i // d for i in vector])
return vector
[docs]def get_struct_from_CIF(cif_fn):
"""
Get structure from CIF file.
:type cif_fn: str
:param cif_fn: CIF file name
:rtype: `schrodinger.structure.Structure`
:return: Structure from CIF file
:raise ValueError: If file doesn't exist or the structure could not be read
"""
if mmutil.feature_flag_is_enabled(mmutil.NATIVE_MMCIF_READER):
try:
return structure.Structure.read(cif_fn)
except OSError:
raise ValueError('Could not convert %s file. Does file exist?' %
cif_fn)
except infrastructure.PyStructureReaderException:
raise ValueError('Bad CIF file: %s.' % cif_fn)
pdb = fileutils.get_basename(cif_fn)
pdb_fn = pdb + '_out.pdb'
cmd = ['obabel', '-icif', cif_fn, '-xo', '-opdb', '-O', pdb_fn]
subprocess.call(cmd)
try:
return structure.Structure.read(pdb_fn)
except StopIteration:
# This happens when file doesn't exist
raise ValueError('Could not convert %s file. Does file exist?' % cif_fn)
[docs]def translate_atoms(cell, origin=None):
"""
Translate all atoms to the first unit cell.
:type cell: schrodinger.structure.Structure
:param cell: the structure to translate, must have chorus
properties defined
:type origin: list or None
:param origin: the list with origin of the cell to which to translate in
fractional coordinates. If none it will be set to [-0.5, -0.5, -0.5]
for system with s_ffio_ct_type set and ParserWrapper.ORIGIN for others.
"""
# get chorus props
chorus_properties = get_chorus_properties(cell)
# sync with PDB props
set_pbc_properties(cell, chorus_properties)
# get origin
if origin is None:
if cell.property.get('s_ffio_ct_type'):
origin = [-0.5, -0.5, -0.5]
else:
origin = ParserWrapper.ORIGIN
# translate
translate_to_cell(cell, origin=origin)
[docs]def translate_molecules(cell,
centroids=False,
maxes=False,
fract_offset=FRACT_OFFSET):
"""
Translate all molecules to the first unit cell.
:type cell: schrodinger.structure.Structure
:param cell: the structure to translate, must have chorus
properties defined
:param bool centroids: If True, translation is using molecular centroid
:param bool maxes: If True, translation is using maximum in x, y, z direction
for each molecule. Cannot be used together with centroids option
:type fract_offset: float
:param fract_offset: The threshold used to compare floating point
fractional coordinate values and in particular those that are on
the cell boundary
"""
# One and only one must be True
assert centroids ^ maxes
xyz = cell.getXYZ(copy=False)
vecs = get_vectors_from_chorus(cell)
for idx, mol in enumerate(cell.molecule):
aindices = mol.getAtomIndices()
if centroids:
cxyz = numpy.array(
transform.get_centroid(cell, atom_list=aindices)[:3])
elif maxes:
# Get max values along x, y, z
cxyz = xyz[numpy.array(aindices) - 1].max(axis=0)
end = translate_to_1st_cell(cxyz, vecs, fract_offset=fract_offset)
move_vec = end - cxyz
# If null vector, don't shift anything
if transform.get_vector_magnitude(move_vec) < 0.01:
continue
# Translate only coordinates of the current molecule
xyz[numpy.array(aindices) - 1] += move_vec
[docs]def get_cell_fast(asu, symm_ops):
"""
Get crystal cell with all applied symmetry operators on the ASU.
:type asu: `structure.Structure`
:param asu: Asymmetric unit
:type symm_ops: list
:param symm_ops: List of symmetry operations (4 x 4 matrices)
:rtype: `structure.Structure`
:return: Crystal cell
"""
vecs = get_vectors_from_chorus(asu)
frac = trans_cart_to_frac_from_vecs(asu.getXYZ(), *vecs)
struct = asu.copy()
struct.deleteAtoms(range(1, struct.atom_total + 1))
struct_temp = asu.copy()
for oper in symm_ops:
oper = numpy.asarray(oper)
frac_new = numpy.dot(frac, oper[:3, :3]) + oper[:3, 3].flat
struct_temp.setXYZ(trans_frac_to_cart_from_vecs(frac_new, *vecs))
struct.extend(struct_temp)
return move_atoms_into_cell(struct)
[docs]def get_pbc_origin(st):
"""
Return the PBC origin in Angstrom from the given structure.
:type st: `structure.Structure`
:param st: the structure
:raise KeyError: if no PBC exists
:rtype: numpy.array
:return: the origin in Angstrom
"""
try:
params = get_lattice_param_properties(st)
except KeyError:
try:
params = get_chorus_properties(st)
except KeyError:
raise
else:
params = get_params_from_chorus(params)
abc_vec = sum(get_lattice_vectors(*params))
pbc_position = st.property.get(PBC_POSITION_KEY)
if pbc_position:
if pbc_position.startswith(ANCHOR_PREFIX):
origin = numpy.array(get_carts_from_anchor_string(pbc_position))
else:
centroid = transform.get_centroid(st)[:3]
origin = centroid - 0.5 * abc_vec
else:
origin = numpy.array(ParserWrapper.ORIGIN)
return origin
[docs]def get_pbc_bonded_atom_pairs(st,
along_vector_idxs=None,
exclusive=False,
only_infinite=False):
"""
Return a dictionary of PBC bonded atom pairs along the specified
directions.
:type st: `structure.Structure`
:param st: the structure
:type along_vector_idxs: list or None
:param along_vector_idxs: the directions to consider,
0 for a-vector, 1 for b-vector, and/or 2 for
c-vector, if None then the c-vector is used by default
:type exclusive: bool
:param exclusive: if True then consider PBC bonds that wrap
only the specified directions, if False then must wrap
the specified directions but can also wrap other directions
:type only_infinite: bool
:param only_infinite: if True then consider only PBC bonds that
are responsible for the system being infinite (meaning that
breaking such bonds makes the infinite system finite)
:rtype: dict
:return: contains atom index pairs for PBC bonds, keys are
direction indices, values are lists of atom index (near, far)
pair tuples where the fractional coordinate along the keyed
direction is near < far
"""
if along_vector_idxs is None:
along_vector_idxs = [2]
atom_idxs = list(range(1, st.atom_total + 1))
pbc = infrastructure.PBC(st)
vecs = pbc.getBoxVectors()
fracs = trans_cart_to_frac_from_vecs(st.getXYZ(), *vecs)
# if speed becomes an issue for the only_infinite option the following
# functions are also applicable: analyze.find_shortest_bond_path,
# analyze.group_by_connectivity, clusterstruct.contract_by_molecule,
# and is_infinite2
all_pairs = {}
for idx in along_vector_idxs:
sorted_atom_idxs = sorted(atom_idxs, key=lambda x: fracs[x - 1, idx])
sorted_atom_idxs.reverse()
pairs = []
for atom_idx in sorted_atom_idxs:
if fracs[atom_idx - 1, idx] <= 0.5:
break
atom = st.atom[atom_idx]
for neighbor in atom.bonded_atoms:
image_vec = pbc.hasNearerImage(st, atom_idx, st, neighbor.index)
bond_vec = numpy.array(neighbor.xyz) - numpy.array(atom.xyz)
if image_vec and image_vec[idx] > 0 and bond_vec[idx] < 0:
if exclusive and transform.get_angle_between_vectors(
-bond_vec, vecs[idx]) > 0.1:
continue
if only_infinite:
tst = st.copy()
tst.deleteBond(atom, neighbor)
if tst.mol_total > st.mol_total:
continue
pairs.append((neighbor.index, atom_idx))
all_pairs[idx] = pairs
return all_pairs
[docs]def delete_pbc_bonds(st,
along_vector_idxs=None,
exclusive=False,
only_infinite=False):
"""
Delete PBC bonds from the given structure along the specified
directions.
:type st: `structure.Structure`
:param st: the structure
:type along_vector_idxs: list or None
:param along_vector_idxs: the directions to consider,
0 for a-vector, 1 for b-vector, and/or 2 for
c-vector, if None then the c-vector is used by default
:type exclusive: bool
:param exclusive: if True then consider PBC bonds that wrap
only the specified directions, if False then must wrap
the specified directions but can also wrap other directions
:type only_infinite: bool
:param only_infinite: if True then consider only PBC bonds that
are responsible for the system being infinite (meaning that
breaking such bonds makes the infinite system finite)
:rtype: dict
:return: contains atom index pairs for PBC bonds, keys are
direction indices, values are lists of atom index (near, far)
pair tuples where the fractional coordinate along the keyed
direction is near < far
"""
all_pairs = get_pbc_bonded_atom_pairs(st,
along_vector_idxs=along_vector_idxs,
exclusive=exclusive,
only_infinite=only_infinite)
for pairs in all_pairs.values():
for pair in pairs:
st.deleteBond(*pair)
return all_pairs
[docs]def has_pbc(st):
"""
Return True if the given structure has a PBC.
:type st: `schrodinger.structure.Structure`
:param st: the structure
:rtype: bool
:return: True if the given structure has a PBC
"""
try:
chorus_properties = get_chorus_properties(st)
except KeyError:
try:
lattice_properties = get_lattice_param_properties(st)
except KeyError:
return False
return True
[docs]def get_pymatgen_structure(st):
"""
Return a pymatgen.core.structure.Structure from the given
schrodinger.structure.Structure.
:type st: `schrodinger.structure.Structure`
:param st: the structure
:raise ValueError: if there is an issue with the input
:rtype: pymatgen.core.structure.Structure
:return: the structure
"""
# This takes approx 2.5 seconds to import
from pymatgen.core.structure import Structure as PMGStructure
if not has_pbc(st):
raise ValueError('No PBC found.')
try:
chorus_to_lower_triangle(st)
chorus_properties = get_chorus_properties(st)
except (ValueError, KeyError):
lattice_properties = get_lattice_param_properties(st)
else:
lattice_properties = get_params_from_chorus(chorus_properties)
vecs = get_lattice_vectors(*lattice_properties)
ans = []
for atom in st.atom:
an = atom.atomic_number
if an < 1:
raise ValueError('Dummy atoms currently not supported.')
ans.append(an)
carts = st.getXYZ()
return PMGStructure(vecs,
ans,
carts,
charge=st.formal_charge,
validate_proximity=False,
to_unit_cell=False,
coords_are_cartesian=True,
site_properties=None)
[docs]def scale_cell(vecs, new_volume):
"""
Return a new Lattice with volume new_volume by performing a scaling of the
lattice vectors so that length proportions and angles are preserved.
:param numpy.array vecs: Original lattice vectors
:param float new_volume: New volume
:rtype: numpy.array
:return: New lattice vectors with the desired volume
"""
# Originally from pymatgen: pymatgen/core/lattice.py :: Lattice::scale
abc = numpy.array(numpy.sqrt(numpy.sum(vecs**2, axis=1)).tolist())
versors = vecs / abc
geo_factor = abs(numpy.dot(numpy.cross(versors[0], versors[1]), versors[2]))
ratios = abc / abc[2]
new_c = (new_volume / (geo_factor * numpy.prod(ratios)))**(1 / 3.0)
return versors * (new_c * ratios)
[docs]def both_bond_types_exist(st):
"""
Raise a ValueError if the given structure does not have at least both
a single regular bond and a single PBC bond, otherwise return sample
atom indices for both.
:type st: `schrodinger.structure.Structure`
:param st: the structure to check
:raise ValueError: if both bond types do not exist
:rtype: tuple
:return: pair tuple containing two pair tuples, one for a regular
bond and one for a PBC bond
"""
pbc_bond_idxs, reg_bond_idxs = (), ()
for bond in st.bond:
if is_pbc_bond(st, bond.atom1, bond.atom2)[0]:
pbc_bond_idxs = (bond.atom1.index, bond.atom2.index)
else:
reg_bond_idxs = (bond.atom1.index, bond.atom2.index)
if pbc_bond_idxs and reg_bond_idxs:
break
else:
raise ValueError('Structure does not have both a regular bond and a '
'PBC bond.')
return tuple((reg_bond_idxs, pbc_bond_idxs))
[docs]def strain_cell(struct, vecs, remap=True):
"""
Strain cell using the vectors and atom coordinates are scaled proportionally
if remap is True.
:param structure.Structure struct: Structure to strain
:param bool remap: whether to strain the atom coordinates in the cell
:param list[3x3] new_vecs: New lattice vectors
"""
vecs = numpy.asarray(vecs)
if not remap:
set_pbc_properties(struct, vecs.flatten())
return
orig_vecs = get_vectors_from_chorus(struct)
frac = trans_cart_to_frac_from_vecs(struct.getXYZ(), *orig_vecs)
cart = trans_frac_to_cart_from_vecs(frac, *vecs)
struct.setXYZ(cart)
set_pbc_properties(struct, vecs.flatten())
[docs]def chorus_to_lower_triangle(struct):
"""
Update chorus properties and with them atom coordinates so that upper
triangle of the lattice vectors are all zeros. This is needed so that
maestro displays structure correctly.
:param structure.Structure struct: Structure to update
:rtype: structure.Structure
:return: Updated structure (not a copy!!)
"""
upper_triangle_indices = [1, 2, 5]
vecs = get_vectors_from_chorus(struct)
if not vecs.flat[upper_triangle_indices].any():
return struct
fracs = trans_cart_to_frac_from_vecs(struct.getXYZ(), *vecs)
params = get_params_from_chorus(vecs.flat)
prim_vecs = numpy.array(get_chorus_from_params(params)).reshape(3, 3)
coords = trans_frac_to_cart_from_vecs(fracs, *prim_vecs)
struct.setXYZ(coords)
set_pbc_properties(struct, prim_vecs.flat)
return struct