diffpy.Structure.SymmetryUtilities

1 ############################################################################## 2 # 3 # Structure by DANSE Diffraction group 4 # Simon J. L. Billinge 5 # (c) 2006 trustees of the Michigan State University. 6 # All rights reserved. 7 # 8 # File coded by: Pavol Juhas 9 # 10 # See AUTHORS.txt for a list of people who contributed. 11 # See LICENSE.txt for license information. 12 # 13 ############################################################################## 14 15 """Symmetry utility functions such as expansion of asymmetric unit, 16 and generation of positional constraints. 17 """ 18 19 # version 20 __id__ = '$Id: SymmetryUtilities.py 2896 2009-03-16 05:45:06Z juhas $' 21 22 import sys 23 import re 24 import numpy 25 26 # Constants: 27 28 # Default tolerance for equality of 2 positions, also 29 # used for identification of special positions. 30 epsilon = 1.0e-5 31 32 # Standard symbols denoting elements of anisotropic thermal 33 # displacement tensor. 34 stdUsymbols = ['U11', 'U22', 'U33', 'U12', 'U13', 'U23'] 35 36 # End of constants 37

38 -def isSpaceGroupLatPar(spacegroup, a, b, c, alpha, beta, gamma):

39 """Check if space group allows passed lattice parameters. 40 41 spacegroup -- instance of SpaceGroup 42 a, b, c, alpha, beta, gamma -- lattice parameters 43 44 Return bool. 45 """ 46 # crystal system rules 47 # ref: Benjamin, W. A., Introduction to crystallography, 48 # New York (1969), p.60 49 def check_triclinic(): 50 return True

51 def check_monoclinic(): 52 rv = (alpha == gamma == 90) or (alpha == beta == 90) 53 return rv 54 def check_orthorhombic(): 55 return (alpha == beta == gamma == 90) 56 def check_tetragonal(): 57 return (a == b and alpha == beta == gamma == 90) 58 def check_trigonal(): 59 rv = (a == b == c and alpha == beta == gamma) or \ 60 (a == b and alpha == beta == 90 and gamma == 120) 61 return rv 62 def check_hexagonal(): 63 return (a == b and alpha == beta == 90 and gamma == 120) 64 def check_cubic(): 65 return (a == b == c and alpha == beta == gamma == 90) 66 crystal_system_rules = { 67 "TRICLINIC" : check_triclinic, 68 "MONOCLINIC" : check_monoclinic, 69 "ORTHORHOMBIC" : check_orthorhombic, 70 "TETRAGONAL" : check_tetragonal, 71 "TRIGONAL" : check_trigonal, 72 "HEXAGONAL" : check_hexagonal, 73 "CUBIC" : check_cubic 74 } 75 rule = crystal_system_rules[spacegroup.crystal_system] 76 return rule() 77 78 # End of isSpaceGroupLatPar 79 80 81 # Constant regular expression used in isconstantFormula(). 82 # isconstantFormula runs faster when regular expression is not 83 # compiled per every single call. 84 85 _rx_constant_formula = re.compile( 86 '[-+]?(\d+(\.\d*)?|\.\d+)([eE][-+]?\d+)??(/[-+]?\d+)?$') 87

88 -def isconstantFormula(s):

89 """Check if formula string is constant. True when argument 90 is a floating point number or a fraction of float with integer. 91 92 s -- formula string 93 94 Return bool. 95 """ 96 res = _rx_constant_formula.match(s.replace(' ', '')) 97 return bool(res)

98 99 # End of isconstantFormula 100 101 102 # Helper class intended for this module only: 103

104 -class _Position2Tuple:

105 """Create callable object that converts fractional coordinates to 106 a tuple of integers with given precision. For presision close to zero 107 it will return a tuples of double. 108 109 Note: Helper class intended for local use only. 110 111 Data members: 112 113 eps -- cutoff for equivalent coordinates. When two coordiantes map to the 114 same tuple, they are closer than eps. 115 """

116 - def __init__(self, eps=epsilon):

117 """Initialize _Position2Tuple 118 119 eps -- cutoff for equivalent coordinates 120 """ 121 # ensure self.eps has exact machine representation 122 self.eps = eps + 1.0 123 self.eps = self.eps - 1.0 124 # no conversions for very small eps 125 if self.eps == 0.0 or 1.0/self.eps > sys.maxint: 126 self.eps = 0.0 127 return

128

129 - def __call__(self, xyz):

130 """Convert array of fractional coordinates to a tuple. 131 132 xyz -- fractional coordinates 133 134 Return a tuple of 3 numbers. 135 """ 136 # no conversion case 137 if self.eps == 0.0: 138 tpl = tuple(xyz % 1.0) 139 return tpl 140 # here we convert to integer 141 tpl = tuple( [int((xi - numpy.floor(xi))/self.eps) for xi in xyz] ) 142 return tpl

143 144 # End of class _Position2Tuple 145

146 -def positionDifference(xyz0, xyz1):

147 """Smallest difference between two coordinates in periodic lattice. 148 149 xyz0, xyz1 -- fractional coordinates 150 151 Return dxyz, a numpy.array dxyz with 0 <= dxyz <= 0.5. 152 """ 153 dxyz = numpy.array(xyz0) - xyz1 154 # map differences to [0,0.5] 155 dxyz = dxyz - numpy.floor(dxyz) 156 mask = (dxyz > 0.5) 157 dxyz[mask] = 1.0 - dxyz[mask] 158 return dxyz

159 160 # End of positionDifference 161

162 -def nearestSiteIndex(sites, xyz):

163 """Index of the nearest site to a specified position. 164 165 sites -- list of coordinates or a 2-dimensional numpy.array 166 xyz -- single position 167 168 Return integer. 169 """ 170 # we use box distance to be consistent with _Position2Tuple conversion 171 dbox = positionDifference(sites, xyz).max(axis=1) 172 nearindex = numpy.argmin(dbox) 173 return nearindex

174 175 # End of nearestSiteIndex 176

177 -def equalPositions(xyz0, xyz1, eps=epsilon):

178 """Equality of two coordinates with optional tolerance. 179 180 xyz0, xyz1 -- fractional coordinates 181 eps -- tolerance for equality of coordinates 182 183 Return bool. 184 """ 185 # we use box distance to be consistent with _Position2Tuple conversion 186 dxyz = positionDifference(xyz0, xyz1) 187 return numpy.all(dxyz <= eps)

188 189 # End of equalPositions 190

191 -def expandPosition(spacegroup, xyz, sgoffset=[0,0,0], eps=epsilon):

192 """Obtain unique equivalent positions and corresponding operations. 193 194 spacegroup -- instance of SpaceGroup 195 xyz -- position to be expanded 196 sgoffset -- offset of space group origin [0, 0, 0] 197 eps -- cutoff for equal positions 198 199 Return a tuple with (list of unique equivalent positions, nested 200 list of SpaceGroups.SymOp instances, site multiplicity). 201 """ 202 sgoffset = numpy.array(sgoffset, dtype=float) 203 pos2tuple = _Position2Tuple(eps) 204 positions = [] 205 site_symops = {} # position tuples with [related symops] 206 for symop in spacegroup.iter_symops(): 207 # operate on coordinates in non-shifted spacegroup 208 pos = symop(xyz + sgoffset) - sgoffset 209 mask = numpy.logical_or(pos < 0.0, pos >= 1.0) 210 pos[mask] -= numpy.floor(pos[mask]) 211 tpl = pos2tuple(pos) 212 if not tpl in site_symops: 213 pos_is_new = True 214 site_symops[tpl] = [] 215 # double check if there is any position nearby 216 if positions: 217 nearpos = positions[nearestSiteIndex(positions, pos)] 218 # is it an equivalent position? 219 if equalPositions(nearpos, pos, eps): 220 # tpl should map to the same list as nearpos 221 site_symops[tpl] = site_symops[ pos2tuple(nearpos) ] 222 pos_is_new = False 223 if pos_is_new: positions.append(pos) 224 # here tpl is inside site_symops 225 site_symops[tpl].append(symop) 226 # pos_symops is nested list of symops associated with each position 227 pos_symops = [ site_symops[pos2tuple(pos)] for pos in positions ] 228 multiplicity = len(positions) 229 return positions, pos_symops, multiplicity

230 231 # End of expandPosition 232

233 -def nullSpace(A):

234 """Null space of matrix A. 235 """ 236 from numpy import linalg 237 u, s, v = linalg.svd(A) 238 # s may have smaller dimension than v 239 vnrows = numpy.shape(v)[0] 240 mask = numpy.ones(vnrows, dtype=bool) 241 mask[s > epsilon] = False 242 null_space = numpy.compress(mask, v, axis=0) 243 return null_space

244 245 # End of nullSpace 246 247

248 -def _findInvariants(symops):

249 """Find a list of symmetry operations which contains identity. 250 251 symops -- nested list of SymOp instances 252 253 Return the list-item in symops which contains identity. 254 Raise ValueError when identity was not found. 255 """ 256 invrnts = None 257 R0 = numpy.identity(3, dtype=float) 258 t0 = numpy.zeros(3, dtype=float) 259 for ops in symops: 260 for op in ops: 261 if numpy.all(op.R == R0) and numpy.all(op.t == t0): 262 invrnts = ops 263 break 264 if invrnts: break 265 if invrnts is None: 266 emsg = "Could not find identity operation." 267 raise ValueError(emsg) 268 return invrnts

269 270 # End of _findInvariants 271 272

273 -class GeneratorSite:

274 """Storage of data related to a generator positions. 275 276 Data members: 277 xyz -- fractional coordinates of generator site 278 Uij -- anisotropic thermal displacement at generator site 279 sgoffset -- offset of space group origin [0, 0, 0] 280 eps -- cutoff for equal positions 281 eqxyz -- list of equivalent positions 282 eqUij -- list of displacement matrices at equivalent positions 283 symops -- nested list of operations per each eqxyz 284 multiplicity -- generator site multiplicity 285 Uisotropy -- bool flag for isotropic thermal factors 286 invariants -- list of invariant operations for generator site 287 null_space -- null space of all possible differences of invariant 288 rotational matrices, this is a base of symmetry 289 allowed shifts. 290 Uspace -- 3D array of independent components of U matrices. 291 pparameters -- list of (xyz symbol, value) pairs 292 Uparameters -- list of (U symbol, value) pairs 293 """ 294 295 Ucomponents = numpy.array([ 296 [[1, 0, 0], [0, 0, 0], [0, 0, 0]], 297 [[0, 0, 0], [0, 1, 0], [0, 0, 0]], 298 [[0, 0, 0], [0, 0, 0], [0, 0, 1]], 299 [[0, 1, 0], [1, 0, 0], [0, 0, 0]], 300 [[0, 0, 1], [0, 0, 0], [1, 0, 0]], 301 [[0, 0, 0], [0, 0, 1], [0, 1, 0]],], dtype=float) 302 idx2Usymbol = { 0 : 'U11', 1 : 'U12', 2 : 'U13', 303 3 : 'U12', 4 : 'U22', 5 : 'U23', 304 6 : 'U13', 7 : 'U23', 8 : 'U33' } 305

306 - def __init__(self, spacegroup, xyz, Uij=numpy.zeros((3,3)), 307 sgoffset=[0,0,0], eps=epsilon):

308 """Initialize GeneratorSite. 309 310 spacegroup -- instance of SpaceGroup 311 xyz -- generating site. When xyz is close to special 312 position self.xyz will be adjusted. 313 Uij -- thermal factors at generator site. Yields self.Uij 314 after adjusting to spacegroup symmetry. 315 sgoffset -- offset of space group origin [0, 0, 0] 316 eps -- cutoff for equal positions 317 """ 318 # just declare the members 319 self.xyz = None 320 self.Uij = Uij 321 self.sgoffset = numpy.array(sgoffset, dtype=float) 322 self.eps = eps 323 self.eqxyz = [] 324 self.eqUij = [] 325 self.symops = None 326 self.multiplicity = None 327 self.Uisotropy = False 328 self.invariants = [] 329 self.null_space = None 330 self.Uspace = None 331 self.pparameters = [] 332 self.Uparameters = [] 333 # fill in the values 334 self.xyz = xyz 335 sites, ops, mult = expandPosition(spacegroup, xyz, sgoffset, eps) 336 invariants = _findInvariants(ops) 337 # shift self.xyz exactly to the special position 338 if mult > 1: 339 xyzdups = numpy.array([op(xyz + self.sgoffset) - self.sgoffset 340 for op in invariants]) 341 dxyz = xyzdups - xyz 342 dxyz = numpy.mean(dxyz - dxyz.round(), axis=0) 343 # recalculate if needed 344 if numpy.any(dxyz != 0.0): 345 self.xyz = xyz + dxyz 346 self.xyz[numpy.fabs(self.xyz) < self.eps] = 0.0 347 sites, ops, mult = expandPosition(spacegroup, 348 self.xyz, self.sgoffset, eps) 349 invariants = _findInvariants(ops) 350 # self.xyz, sites, ops are all adjusted here 351 self.eqxyz = sites 352 self.symops = ops 353 self.multiplicity = mult 354 self.invariants = invariants 355 self._findNullSpace() 356 self._findPosParameters() 357 self._findUSpace() 358 self._findUParameters() 359 self._findeqUij() 360 return

361

362 - def signedRatStr(self, x):

363 """Convert floating point number to signed rational representation. 364 Possible fractional are multiples of 1/3, 1/6, 1/7, 1/9, if these 365 are not close, return "%+g" format. 366 367 Return string. 368 """ 369 s = str(x) 370 if len(s) < 6: return "%+g" % x 371 den = numpy.array([3.0, 6.0, 7.0, 9.0]) 372 nom = x * den 373 idx = numpy.where(numpy.fabs(nom - nom.round()) < self.eps)[0] 374 if idx.size == 0: return "%+g" % x 375 # here we have fraction 376 return "%+.0f/%.0f" % (nom[idx[0]], den[idx[0]])

377

378 - def _findNullSpace(self):

379 """Calculate self.null_space from self.invariants. 380 Try to represent self.null_space using small integers. 381 """ 382 R0 = self.invariants[0].R 383 Rdiff = [ (symop.R - R0) for symop in self.invariants ] 384 Rdiff = numpy.concatenate(Rdiff, 0) 385 self.null_space = nullSpace(Rdiff) 386 if self.null_space.size == 0: return 387 # reverse sort rows of null_space rows by absolute value 388 key = tuple(numpy.fabs(numpy.transpose(self.null_space))[::-1]) 389 order = numpy.lexsort(key) 390 self.null_space = self.null_space[order[::-1]] 391 # rationalize by the smallest element larger than cutoff 392 cutoff = 1.0/32 393 for i in range(len(self.null_space)): 394 row = self.null_space[i] 395 small = numpy.fabs(row[numpy.fabs(row) > cutoff]).min() 396 signedsmall = row[numpy.fabs(row) == small][0] 397 self.null_space[i] = self.null_space[i] / signedsmall 398 return

399

400 - def _findPosParameters(self):

401 """Find pparameters and their values for expressing self.xyz. 402 """ 403 usedsymbol = {} 404 # parameter values depend on offset of self.xyz 405 txyz = self.xyz 406 # define txyz such that most of its elements are zero 407 for nvec in self.null_space: 408 idx = numpy.where(numpy.fabs(nvec) >= epsilon)[0][0] 409 varvalue = txyz[idx]/nvec[idx] 410 txyz = txyz - varvalue*nvec 411 # determine standard parameter name 412 vname = [s for s in "xyz"[idx:] if not s in usedsymbol][0] 413 self.pparameters.append( (vname, varvalue) ) 414 usedsymbol[vname] = True 415 return

416

417 - def _findUSpace(self):

418 """Find independent U components with respect to invariant 419 rotations. 420 """ 421 self.Uspace = [] 422 independent = dict.fromkeys(range(6)) 423 for idx in range(6): 424 if idx not in independent: continue 425 # new U space matrix 426 Uc = self.Ucomponents[idx] 427 Usp = numpy.zeros((3, 3), dtype=float) 428 for op in self.invariants: 429 R = op.R 430 Rt = R.transpose() 431 mxrot = lambda M : numpy.dot(R, numpy.dot(M, Rt)) 432 opUsp = numpy.zeros((3, 3), dtype=float) 433 uceqnext = Uc 434 while True: 435 opUsp += uceqnext 436 uceqnext = mxrot(uceqnext) 437 # does uceqnext contain -Uc? 438 # if yes, the component Uc is not in the Uspace 439 if numpy.all(uceqnext*Uc == -Uc): 440 opUsp[:] = 0.0 441 Usp[:] = 0.0 442 break 443 # all Uc rotations were visited when we are back to Uc 444 if numpy.all(uceqnext == Uc): 445 break 446 found_superspace = (numpy.all(opUsp[Usp != 0]) and 447 numpy.any(opUsp[Usp == 0])) 448 if found_superspace: 449 Usp = opUsp 450 # Usp is all zero here only when component Uc is not 451 # in the Uspace. Do not check the remaining invariants. 452 if numpy.all(Usp == 0.0): break 453 # Find which U components are contained in Usp 454 Usp_contains = [i for i in independent 455 if numpy.any(Usp * self.Ucomponents[i])] 456 # Usp must be all zero if it does not contain any component 457 if not Usp_contains: 458 assert numpy.all(Usp == 0.0) 459 independent.pop(idx) 460 continue 461 # U components contained in Usp are not independent anymore. 462 for i in Usp_contains: del independent[i] 463 # orthogonalize Usp with respect to preceding Uspace components 464 Uspflat = Usp.flatten() 465 for Usp0 in self.Uspace: 466 Usp0flat = Usp0.flatten() 467 projection = numpy.dot(Uspflat, Usp0flat) 468 if projection != 0: 469 projection /= numpy.dot(Usp0flat, Usp0flat) 470 Uspflat -= projection * Usp0flat 471 # normalize Uspflat by its maximum component 472 maxidx = numpy.argmax(numpy.fabs(Uspflat)) 473 Uspflat /= Uspflat[maxidx] 474 Usp = Uspflat.reshape(3, 3) 475 self.Uspace.append(Usp) 476 # finally convert this to 3D array 477 self.Uspace = numpy.array(self.Uspace) 478 self.Uisotropy = (len(self.Uspace) == 1) 479 return

480

481 - def _findUParameters(self):

482 """Find Uparameters and their values for expressing self.Uij. 483 """ 484 # permute indices as 00 11 22 01 02 12 10 20 21 485 diagorder = numpy.array((0, 4, 8, 1, 2, 5, 3, 6, 7)) 486 Uijflat = self.Uij.flatten() 487 for Usp in self.Uspace: 488 Uspflat = Usp.flatten() 489 Uspnorm2 = numpy.dot(Uspflat, Uspflat) 490 permidx = numpy.where(Uspflat[diagorder])[0][0] 491 idx = diagorder[permidx] 492 vname = self.idx2Usymbol[idx] 493 varvalue = numpy.dot(Uijflat, Uspflat) / Uspnorm2 494 self.Uparameters.append( (vname, varvalue) ) 495 return

496

497 - def _findeqUij(self):

498 """Adjust self.Uij and self.eqUij to be consistent with spacegroup 499 """ 500 self.Uij = numpy.zeros((3,3), dtype=float) 501 for i in range(len(self.Uparameters)): 502 Usp = self.Uspace[i] 503 varvalue = self.Uparameters[i][1] 504 self.Uij += varvalue*Usp 505 # now determine eqUij 506 for ops in self.symops: 507 # take first rotation matrix 508 R = ops[0].R 509 Rt = R.transpose() 510 self.eqUij.append( numpy.dot(R, numpy.dot(self.Uij, Rt)) ) 511 return

512

513 - def positionFormula(self, pos, xyzsymbols=("x","y","z")):

514 """Formula of equivalent position with respect to generator site 515 516 pos -- fractional coordinates of possibly equivalent site 517 xyzsymbols -- symbols for parametrized coordinates 518 519 Return position formulas in a dictionary with keys equal ("x", "y", "z") 520 or an empty dictionary when pos is not equivalent to generator. 521 Formulas are formatted as "[[-][%g*]{x|y|z}] [{+|-}%g]", for example 522 "-x", "z +0.5", "0.25". 523 """ 524 # find pos in eqxyz 525 idx = nearestSiteIndex(self.eqxyz, pos) 526 eqpos = self.eqxyz[idx] 527 if not equalPositions(eqpos, pos, self.eps): return {} 528 # any rotation matrix should do fine 529 R = self.symops[idx][0].R 530 nsrotated = numpy.dot(self.null_space, numpy.transpose(R)) 531 # build formulas using eqpos 532 # find offset 533 teqpos = numpy.array(eqpos) 534 for nvec, (vname, varvalue) in zip(nsrotated, self.pparameters): 535 teqpos -= nvec * varvalue 536 # map varnames to xyzsymbols 537 name2sym = dict( zip(("x", "y", "z"), xyzsymbols) ) 538 xyzformula = 3*[""] 539 for nvec, (vname, ignore) in zip(nsrotated, self.pparameters): 540 for i in range(3): 541 if abs(nvec[i]) < epsilon: continue 542 xyzformula[i] += "%s*%s " % \ 543 (self.signedRatStr(nvec[i]), name2sym[vname]) 544 # add constant offset teqpos to all formulas 545 for i in range(3): 546 if xyzformula[i] and abs(teqpos[i]) < epsilon: continue 547 xyzformula[i] += self.signedRatStr(teqpos[i]) 548 # reduce unnecessary +1* and -1* 549 xyzformula = [ re.sub('^[+]1[*]|(?<=[+-])1[*]', '', f).strip() 550 for f in xyzformula ] 551 return dict( zip(("x","y","z"), xyzformula) )

552

553 - def UFormula(self, pos, Usymbols=stdUsymbols):

554 """List of atom displacement formulas with custom parameter symbols. 555 556 pos -- fractional coordinates of possibly equivalent site 557 Usymbols -- 6 symbols for possible U matrix parameters 558 559 Return U element formulas in a dictionary where keys are from 560 ('U11','U22','U33','U12','U13','U23') or empty dictionary when 561 pos is not equivalent to generator. 562 """ 563 # find pos in eqxyz 564 idx = nearestSiteIndex(self.eqxyz, pos) 565 eqpos = self.eqxyz[idx] 566 if not equalPositions(eqpos, pos, self.eps): return {} 567 # any rotation matrix should do fine 568 R = self.symops[idx][0].R 569 Rt = R.transpose() 570 Usrotated = [numpy.dot(R, numpy.dot(Us, Rt)) for Us in self.Uspace] 571 Uformula = dict.fromkeys(stdUsymbols, '0') 572 name2sym = dict(zip(stdUsymbols, Usymbols)) 573 for Usr, (vname, ignore) in zip(Usrotated, self.Uparameters): 574 Usrflat = Usr.flatten() 575 for i in numpy.where(Usrflat)[0]: 576 f = '%g*%s' % (Usrflat[i], name2sym[vname]) 577 f = re.sub('^[+]?1[*]|(?<=[+-])1[*]', '', f).strip() 578 smbl = self.idx2Usymbol[i] 579 Uformula[smbl] = f 580 return Uformula

581

582 - def eqIndex(self, pos):

583 """Index of the nearest generator equivalent site 584 585 pos -- fractional coordinates 586 587 Return integer. 588 """ 589 return nearestSiteIndex(self.eqxyz, pos)

590 591 # End of class GeneratorSite 592

593 -class ExpandAsymmetricUnit:

594 """Expand asymmetric unit and anisotropic thermal displacement 595 596 Data members: 597 spacegroup -- instance of SpaceGroup 598 corepos -- list of positions in asymmetric unit, 599 it may contain duplicates 600 coreUijs -- thermal factors for corepos (defaults to zeros) 601 sgoffset -- optional offset of space group origin [0, 0, 0] 602 eps -- cutoff for equivalent positions 603 Calculated data members: 604 multiplicity -- multiplicity of each site in corepos 605 Uisotropy -- bool flags for isotropic sites in corepos 606 expandedpos -- list of equivalent positions per each site in corepos 607 expandedUijs -- list of thermal factors per each site in corepos 608 """ 609 610 # By design Atom instances are not accepted as arguments to keep 611 # number of required imports low.

612 - def __init__(self, spacegroup, corepos, coreUijs=None, 613 sgoffset=[0,0,0], eps=epsilon):

614 """Initialize and calculate instance of ExpandAsymmetricUnit 615 616 spacegroup -- instance of SpaceGroup 617 corepos -- list of positions in asymmetric unit, 618 it may contain duplicates 619 coreUijs -- thermal factors for corepos (defaults to zeros) 620 sgoffset -- offset of space group origin [0, 0, 0] 621 eps -- cutoff for duplicate positions 622 """ 623 # declare data members 624 self.spacegroup = spacegroup 625 self.corepos = corepos 626 self.coreUijs = None 627 self.sgoffset = numpy.array(sgoffset) 628 self.eps = eps 629 self.multiplicity = [] 630 self.Uisotropy = [] 631 self.expandedpos = [] 632 self.expandedUijs = [] 633 # obtain their values 634 corelen = len(self.corepos) 635 if coreUijs: 636 self.coreUijs = coreUijs 637 else: 638 self.coreUijs = numpy.zeros((corelen,3,3), dtype=float) 639 for cpos, cUij in zip(self.corepos, self.coreUijs): 640 gen = GeneratorSite(self.spacegroup, cpos, cUij, 641 self.sgoffset, self.eps) 642 self.multiplicity.append(gen.multiplicity) 643 self.Uisotropy.append(gen.Uisotropy) 644 self.expandedpos.append(gen.eqxyz) 645 self.expandedUijs.append(gen.eqUij) 646 return

647 648 # End of ExpandAsymmetricUnit 649 650 651 # Helper function for SymmetryConstraints class. It may be useful 652 # elsewhere therefore its name does not start with underscore. 653

654 -def pruneFormulaDictionary(eqdict):

655 """Remove constant items from formula dictionary. 656 657 eqdict -- formula dictionary which maps standard variable symbols 658 ("x", "U11") to string formulas ("0", "-x3", "z7 +0.5") 659 660 Return pruned formula dictionary. 661 """ 662 pruned = {} 663 for smb, eq in eqdict.iteritems(): 664 if not isconstantFormula(eq): pruned[smb] = eq 665 return pruned

666 667 # End of pruneFormulaDictionary 668

669 -class SymmetryConstraints:

670 """Generate symmetry constraints for specified positions 671 672 Data members: 673 spacegroup -- instance of SpaceGroup 674 positions -- all positions to be constrained 675 Uijs -- thermal factors for all positions (defaults to zeros) 676 sgoffset -- optional offset of space group origin [0, 0, 0] 677 eps -- cutoff for equivalent positions 678 Calculated data members: 679 corepos -- list of of positions in the asymmetric unit 680 coremap -- dictionary mapping indices of asymmetric core positions 681 to indices of all symmetry related positions 682 poseqns -- list of coordinate formula dictionaries per each site. 683 Formula dictionary keys are from ("x", "y", "z") and 684 the values are formatted as [[-]{x|y|z}%i] [{+|-}%g], 685 for example: "x0", "-x3", "z7 +0.5", "0.25". 686 pospars -- list of (xyz symbol, value) pairs 687 Ueqns -- list of anisotropic atomic displacement formula 688 dictionaries per each position. Formula dictionary 689 keys are from ('U11','U22','U33','U12','U13','U23') 690 and the values are formatted as {[%g*][Uij%i]|0}, 691 for example: "U110", "0.5*U2213", "0" 692 Upars -- list of (U symbol, value) pairs 693 Uisotropy -- list of bool flags for isotropic thermal displacements 694 """ 695

696 - def __init__(self, spacegroup, positions, Uijs=None, 697 sgoffset=[0, 0, 0], eps=epsilon):

698 """Initialize and calculate SymmetryConstraints. 699 700 spacegroup -- instance of SpaceGroup 701 positions -- list of all positions to be constrained 702 Uijs -- list of U matrices for all constrained positions 703 sgoffset -- optional offset of space group origin [0, 0, 0] 704 eps -- cutoff for equivalent positions 705 """ 706 # fill in data members 707 self.spacegroup = spacegroup 708 self.positions = None 709 self.Uijs = None 710 self.sgoffset = numpy.array(sgoffset) 711 self.eps = eps 712 self.corepos = [] 713 self.coremap = {} 714 self.poseqns = None 715 self.pospars = [] 716 self.Ueqns = None 717 self.Upars = [] 718 self.Uisotropy = None 719 # handle single position: 720 import types 721 # handle list of lists returned by ExpandAsymmetricUnit 722 if len(positions) and type(positions[0]) is types.ListType: 723 # concatenate lists before converting to Nx3 array 724 flatpos = sum(positions, []) 725 flatpos = numpy.array(flatpos, dtype=float).flatten() 726 self.positions = flatpos.reshape((flatpos.size/3, 3)) 727 # otherwise convert to array 728 else: 729 flatpos = numpy.array(positions, dtype=float).flatten() 730 self.positions = flatpos.reshape((flatpos.size/3, 3)) 731 # here self.positions should be a 2D numpy array 732 numpos = len(self.positions) 733 # adjust Uijs if not specified 734 if Uijs is not None: 735 self.Uijs = numpy.array(Uijs) 736 else: 737 self.Uijs = numpy.zeros((numpos, 3, 3), dtype=float) 738 self.poseqns = numpos*[None] 739 self.Ueqns = numpos*[None] 740 self.Uisotropy = numpos*[False] 741 # all members should be initialized here 742 self._findConstraints() 743 return

744

745 - def _findConstraints(self):

746 """Find constraints for positions and anisotropic displacements Uij 747 """ 748 numpos = len(self.positions) 749 # canonical xyzsymbols and Usymbols 750 xyzsymbols = [ smbl+str(i) for i in range(numpos) for smbl in "xyz" ] 751 Usymbols = [smbl+str(i) for i in range(numpos) for smbl in stdUsymbols] 752 independent = dict.fromkeys(range(numpos)) 753 for genidx in range(numpos): 754 if not genidx in independent: continue 755 # it is a generator 756 self.coremap[genidx] = [] 757 genpos = self.positions[genidx] 758 genUij = self.Uijs[genidx] 759 gen = GeneratorSite(self.spacegroup, genpos, genUij, 760 self.sgoffset, self.eps) 761 # append new pparameters if there are any 762 gxyzsymbols = xyzsymbols[3*genidx : 3*(genidx+1)] 763 for k, v in gen.pparameters: 764 smbl = gxyzsymbols["xyz".index(k)] 765 self.pospars.append( (smbl, v) ) 766 gUsymbols = Usymbols[6*genidx : 6*(genidx+1)] 767 for k, v in gen.Uparameters: 768 smbl = gUsymbols[stdUsymbols.index(k)] 769 self.Upars.append( (smbl, v) ) 770 # search for equivalents inside indies 771 indies = independent.keys() 772 indies.sort() 773 for indidx in indies: 774 indpos = self.positions[indidx] 775 formula = gen.positionFormula(indpos, gxyzsymbols) 776 # formula is empty when indidx is independent 777 if not formula: continue 778 # indidx is dependent here 779 del independent[indidx] 780 self.coremap[genidx].append(indidx) 781 self.poseqns[indidx] = formula 782 self.Ueqns[indidx] = gen.UFormula(indpos, gUsymbols) 783 # make sure positions and Uijs are consistent with spacegroup 784 eqidx = gen.eqIndex(indpos) 785 dxyz = gen.eqxyz[eqidx] - indpos 786 self.positions[indidx] += dxyz - dxyz.round() 787 self.Uijs[indidx] = gen.eqUij[eqidx] 788 # all done here 789 coreidx = self.coremap.keys() 790 coreidx.sort() 791 self.corepos = [self.positions[i] for i in coreidx] 792 return

793

794 - def posparSymbols(self):

795 """Return list of standard position parameter symbols. 796 """ 797 return [n for n, v in self.pospars]

798

799 - def posparValues(self):

800 """Return list of position parameters values. 801 """ 802 return [v for n, v in self.pospars]

803

804 - def positionFormulas(self, xyzsymbols=None):

805 """List of position formulas with custom parameter symbols. 806 807 xyzsymbols -- list of custom symbols used in formula strings 808 809 Return list of coordinate formulas dictionaries. Formulas dictionary 810 keys are from ("x", "y", "z") and the values are formatted as 811 [[-]{symbol}] [{+|-}%g], for example: "x0", "-sym", "@7 +0.5", "0.25". 812 """ 813 if not xyzsymbols: return list(self.poseqns) 814 # check xyzsymbols 815 if len(xyzsymbols) < len(self.pospars): 816 emsg = ("Not enough symbols for %i position parameters" % 817 len(self.pospars)) 818 raise SymmetryError(emsg) 819 # build translation dictionary 820 trsmbl = dict(zip(self.posparSymbols(), xyzsymbols)) 821 def translatesymbol(matchobj): 822 return trsmbl[matchobj.group(0)]

823 pat = re.compile(r'\b[xyz]\d+') 824 rv = [] 825 for eqns in self.poseqns: 826 treqns = {} 827 for smbl, eq in eqns.iteritems(): 828 treqns[smbl] = re.sub(pat, translatesymbol, eq) 829 rv.append(treqns) 830 return rv

831

832 - def positionFormulasPruned(self, xyzsymbols=None):

833 """List of position formula dictionaries with constant items removed. 834 See also positionFormulas(). 835 836 xyzsymbols -- list of custom symbols used in formula strings 837 838 Return list of coordinate formula dictionaries. 839 """ 840 rv = [ pruneFormulaDictionary(eqns) 841 for eqns in self.positionFormulas(xyzsymbols) ] 842 return rv

843

844 - def UparSymbols(self):

845 """Return list of standard atom displacement parameter symbols. 846 """ 847 return [n for n, v in self.Upars]

848

849 - def UparValues(self):

850 """Return list of atom displacement parameters values. 851 """ 852 return [v for n, v in self.Upars]

853

854 - def UFormulas(self, Usymbols=None):

855 """List of atom displacement formulas with custom parameter symbols. 856 857 Usymbols -- list of custom symbols used in formula strings 858 859 Return list of atom displacement formula dictionaries per each site. 860 Formula dictionary keys are from ('U11','U22','U33','U12','U13','U23') 861 and the values are formatted as {[%g*][Usymbol]|0}, for example: 862 "U11", "0.5*@37", "0". 863 """ 864 if not Usymbols: return list(self.Ueqns) 865 # check Usymbols 866 if len(Usymbols) < len(self.Upars): 867 emsg = "Not enough symbols for %i U parameters" % len(self.Upars) 868 raise SymmetryError(emsg) 869 # build translation dictionary 870 trsmbl = dict(zip(self.UparSymbols(), Usymbols)) 871 def translatesymbol(matchobj): 872 return trsmbl[matchobj.group(0)]

873 pat = re.compile(r'\bU\d\d\d+') 874 rv = [] 875 for eqns in self.Ueqns: 876 treqns = {} 877 for smbl, eq in eqns.iteritems(): 878 treqns[smbl] = re.sub(pat, translatesymbol, eq) 879 rv.append(treqns) 880 return rv 881

882 - def UFormulasPruned(self, Usymbols=None):

883 """List of atom displacement formula dictionaries with constant items 884 removed. See also UFormulas(). 885 886 Usymbols -- list of custom symbols used in formula strings 887 888 Return list of atom displacement formulas in tuples of 889 (U11, U22, U33, U12, U13, U23). 890 """ 891 rv = [ pruneFormulaDictionary(eqns) 892 for eqns in self.UFormulas(Usymbols) ] 893 return rv

894 895 # End of SymmetryConstraints 896 897 # basic demonstration 898 if __name__ == "__main__": 899 from diffpy.Structure.SpaceGroups import sg100 900 site = [.125, .625, .13] 901 Uij = [[1,2,3],[2,4,5],[3,5,6]] 902 g = GeneratorSite(sg100, site, Uij=Uij) 903 fm100 = g.positionFormula(site) 904 print "g = GeneratorSite(sg100, %r)" % site 905 print "g.positionFormula(%r) = %s" % (site, fm100) 906 print "g.pparameters =", g.pparameters 907 print "g.Uparameters =", g.Uparameters 908 print "g.UFormula(%r) =" % site, g.UFormula(site) 909 910 # End of file 911

Source Code for Module diffpy.Structure.SymmetryUtilities