| 1 | # -*- coding: utf-8 -*- |
|---|
| 2 | ########### SVN repository information ################### |
|---|
| 3 | # $Date: 2018-05-24 17:19:27 +0000 (Thu, 24 May 2018) $ |
|---|
| 4 | # $Author: vondreele $ |
|---|
| 5 | # $Revision: 3401 $ |
|---|
| 6 | # $URL: trunk/imports/G2phase_CIF.py $ |
|---|
| 7 | # $Id: G2phase_CIF.py 3401 2018-05-24 17:19:27Z vondreele $ |
|---|
| 8 | ########### SVN repository information ################### |
|---|
| 9 | ''' |
|---|
| 10 | *Module G2phase_CIF: Coordinates from CIF* |
|---|
| 11 | ------------------------------------------ |
|---|
| 12 | |
|---|
| 13 | Parses a CIF using PyCifRW from James Hester and pulls out the |
|---|
| 14 | structural information. |
|---|
| 15 | |
|---|
| 16 | If a CIF generated by ISODISTORT is encountered, extra information is |
|---|
| 17 | added to the phase entry and constraints are generated. |
|---|
| 18 | |
|---|
| 19 | ''' |
|---|
| 20 | # Routines to import Phase information from CIF files |
|---|
| 21 | from __future__ import division, print_function |
|---|
| 22 | import sys |
|---|
| 23 | import random as ran |
|---|
| 24 | import numpy as np |
|---|
| 25 | import re |
|---|
| 26 | import GSASIIIO as G2IO |
|---|
| 27 | import GSASIIobj as G2obj |
|---|
| 28 | import GSASIIspc as G2spc |
|---|
| 29 | import GSASIIElem as G2elem |
|---|
| 30 | import GSASIIlattice as G2lat |
|---|
| 31 | import GSASIIpy3 as G2p3 |
|---|
| 32 | import GSASIIpath |
|---|
| 33 | GSASIIpath.SetVersionNumber("$Revision: 3401 $") |
|---|
| 34 | import CifFile as cif # PyCifRW from James Hester |
|---|
| 35 | |
|---|
| 36 | class CIFPhaseReader(G2obj.ImportPhase): |
|---|
| 37 | 'Implements a phase importer from a possibly multi-block CIF file' |
|---|
| 38 | def __init__(self): |
|---|
| 39 | super(self.__class__,self).__init__( # fancy way to say ImportPhase.__init__ |
|---|
| 40 | extensionlist=('.CIF','.cif','.mcif'), |
|---|
| 41 | strictExtension=False, |
|---|
| 42 | formatName = 'CIF', |
|---|
| 43 | longFormatName = 'Crystallographic Information File import' |
|---|
| 44 | ) |
|---|
| 45 | |
|---|
| 46 | def ContentsValidator(self, filename): |
|---|
| 47 | fp = open(filename,'r') |
|---|
| 48 | return self.CIFValidator(fp) |
|---|
| 49 | fp.close() |
|---|
| 50 | |
|---|
| 51 | def Reader(self,filename, ParentFrame=None, usedRanIdList=[], **unused): |
|---|
| 52 | self.isodistort_warnings = '' |
|---|
| 53 | self.Phase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 54 | # make sure the ranId is really unique! |
|---|
| 55 | while self.Phase['ranId'] in usedRanIdList: |
|---|
| 56 | self.Phase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 57 | self.MPhase = G2obj.SetNewPhase(Name='new phase') # create a new empty phase dict |
|---|
| 58 | # make sure the ranId is really unique! |
|---|
| 59 | while self.MPhase['ranId'] in usedRanIdList: |
|---|
| 60 | self.MPhase['ranId'] = ran.randint(0,sys.maxsize) |
|---|
| 61 | returnstat = False |
|---|
| 62 | cellitems = ( |
|---|
| 63 | '_cell_length_a','_cell_length_b','_cell_length_c', |
|---|
| 64 | '_cell_angle_alpha','_cell_angle_beta','_cell_angle_gamma',) |
|---|
| 65 | # cellwaveitems = ( |
|---|
| 66 | # '_cell_wave_vector_seq_id', |
|---|
| 67 | # '_cell_wave_vector_x','_cell_wave_vector_y','_cell_wave_vector_z') |
|---|
| 68 | reqitems = ( |
|---|
| 69 | '_atom_site_fract_x', |
|---|
| 70 | '_atom_site_fract_y', |
|---|
| 71 | '_atom_site_fract_z', |
|---|
| 72 | ) |
|---|
| 73 | phasenamefields = ( |
|---|
| 74 | '_chemical_name_common', |
|---|
| 75 | '_pd_phase_name', |
|---|
| 76 | '_chemical_formula_sum' |
|---|
| 77 | ) |
|---|
| 78 | try: |
|---|
| 79 | cf = G2obj.ReadCIF(filename) |
|---|
| 80 | except cif.StarError as msg: |
|---|
| 81 | msg = 'Unreadable cif file\n'+str(msg) |
|---|
| 82 | self.errors = msg |
|---|
| 83 | self.warnings += msg |
|---|
| 84 | return False |
|---|
| 85 | # scan blocks for structural info |
|---|
| 86 | self.errors = 'Error during scan of blocks for datasets' |
|---|
| 87 | str_blklist = [] |
|---|
| 88 | for blk in cf.keys(): |
|---|
| 89 | for r in reqitems+cellitems: |
|---|
| 90 | if r not in cf[blk].keys(): |
|---|
| 91 | break |
|---|
| 92 | else: |
|---|
| 93 | str_blklist.append(blk) |
|---|
| 94 | if not str_blklist: |
|---|
| 95 | selblk = None # no block to choose |
|---|
| 96 | elif len(str_blklist) == 1: # only one choice |
|---|
| 97 | selblk = 0 |
|---|
| 98 | else: # choose from options |
|---|
| 99 | choice = [] |
|---|
| 100 | for blknm in str_blklist: |
|---|
| 101 | choice.append('') |
|---|
| 102 | # accumumlate some info about this phase |
|---|
| 103 | choice[-1] += blknm + ': ' |
|---|
| 104 | for i in phasenamefields: # get a name for the phase |
|---|
| 105 | name = cf[blknm].get(i,'phase name').strip() |
|---|
| 106 | if name is None or name == '?' or name == '.': |
|---|
| 107 | continue |
|---|
| 108 | else: |
|---|
| 109 | choice[-1] += name.strip()[:20] + ', ' |
|---|
| 110 | break |
|---|
| 111 | na = len(cf[blknm].get("_atom_site_fract_x")) |
|---|
| 112 | if na == 1: |
|---|
| 113 | choice[-1] += '1 atom' |
|---|
| 114 | else: |
|---|
| 115 | choice[-1] += ('%d' % na) + ' atoms' |
|---|
| 116 | choice[-1] += ', cell: ' |
|---|
| 117 | fmt = "%.2f," |
|---|
| 118 | for i,key in enumerate(cellitems): |
|---|
| 119 | if i == 3: fmt = "%.f," |
|---|
| 120 | if i == 5: fmt = "%.f" |
|---|
| 121 | choice[-1] += fmt % cif.get_number_with_esd( |
|---|
| 122 | cf[blknm].get(key))[0] |
|---|
| 123 | sg = cf[blknm].get("_symmetry_space_group_name_H-M",'') |
|---|
| 124 | if not sg: sg = cf[blknm].get("_space_group_name_H-M_alt",'') |
|---|
| 125 | if not sg: sg = cf[blknm].get("_space_group_ssg_name",'') |
|---|
| 126 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 127 | if not sg: sg = cf[blknm].get("_space_group.magn_ssg_name",'') |
|---|
| 128 | #how about checking for super/magnetic ones as well? - reject 'X'? |
|---|
| 129 | sg = sg.replace('_','') |
|---|
| 130 | if sg: choice[-1] += ', (' + sg.strip() + ')' |
|---|
| 131 | selblk = G2IO.PhaseSelector(choice,ParentFrame=ParentFrame, |
|---|
| 132 | title= 'Select a phase from one the CIF data_ blocks below',size=(600,100)) |
|---|
| 133 | self.errors = 'Error during reading of selected block' |
|---|
| 134 | #process selected phase |
|---|
| 135 | if selblk is None: |
|---|
| 136 | returnstat = False # no block selected or available |
|---|
| 137 | else: #do space group symbol & phase type first |
|---|
| 138 | blknm = str_blklist[selblk] |
|---|
| 139 | blk = cf[str_blklist[selblk]] |
|---|
| 140 | E = True |
|---|
| 141 | Super = False |
|---|
| 142 | magnetic = False |
|---|
| 143 | moddim = int(blk.get("_cell_modulation_dimension",'0')) |
|---|
| 144 | if moddim: #incommensurate |
|---|
| 145 | if moddim > 1: |
|---|
| 146 | msg = 'more than 3+1 super symmetry is not allowed in GSAS-II' |
|---|
| 147 | self.errors = msg |
|---|
| 148 | self.warnings += '\n'+msg |
|---|
| 149 | return False |
|---|
| 150 | if blk.get('_cell_subsystems_number'): |
|---|
| 151 | msg = 'Composite super structures not allowed in GSAS-II' |
|---|
| 152 | self.errors = msg |
|---|
| 153 | self.warnings += '\n'+msg |
|---|
| 154 | return False |
|---|
| 155 | sspgrp = blk.get("_space_group_ssg_name",'') |
|---|
| 156 | if not sspgrp: #might be incommensurate magnetic |
|---|
| 157 | MSSpGrp = blk.get("_space_group.magn_ssg_name_BNS",'') |
|---|
| 158 | if not MSSpGrp: |
|---|
| 159 | MSSpGrp = blk.get("_space_group.magn_ssg_name",'') |
|---|
| 160 | if not MSSpGrp: |
|---|
| 161 | msg = 'No incommensurate space group name was found in the CIF.' |
|---|
| 162 | self.errors = msg |
|---|
| 163 | self.warnings += '\n'+msg |
|---|
| 164 | return False |
|---|
| 165 | if 'X' in MSSpGrp: |
|---|
| 166 | msg = 'Ad hoc incommensurate magnetic space group '+MSSpGrp+' is not allowed in GSAS-II' |
|---|
| 167 | self.warnings += '\n'+msg |
|---|
| 168 | self.errors = msg |
|---|
| 169 | return False |
|---|
| 170 | magnetic = True |
|---|
| 171 | if 'X' in sspgrp: |
|---|
| 172 | msg = 'Ad hoc incommensurate space group '+sspgrp+' is not allowed in GSAS-II' |
|---|
| 173 | self.warnings += '\n'+msg |
|---|
| 174 | self.errors = msg |
|---|
| 175 | return False |
|---|
| 176 | Super = True |
|---|
| 177 | if magnetic: |
|---|
| 178 | sspgrp = MSSpGrp.split('(') |
|---|
| 179 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 180 | SpGrp = G2spc.StandardizeSpcName(sspgrp[0]) |
|---|
| 181 | if "1'" in SpGrp: sspgrp[1] = sspgrp[1][:-1] #take off extra 's'; gets put back later |
|---|
| 182 | MSpGrp = sspgrp[0] |
|---|
| 183 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 184 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 185 | else: |
|---|
| 186 | sspgrp = sspgrp.split('(') |
|---|
| 187 | sspgrp[1] = "("+sspgrp[1] |
|---|
| 188 | SpGrp = sspgrp[0] |
|---|
| 189 | SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 190 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 191 | if not SpGrp: |
|---|
| 192 | print (sspgrp) |
|---|
| 193 | self.warnings += 'Space group name '+sspgrp[0]+sspgrp[1]+' not recognized by GSAS-II' |
|---|
| 194 | return False |
|---|
| 195 | SuperSg = sspgrp[1].replace('\\','') |
|---|
| 196 | SuperVec = [[0,0,.1],False,4] |
|---|
| 197 | else: #not incommensurate |
|---|
| 198 | SpGrp = blk.get("_symmetry_space_group_name_H-M",'') |
|---|
| 199 | if not SpGrp: |
|---|
| 200 | SpGrp = blk.get("_space_group_name_H-M_alt",'') |
|---|
| 201 | if not SpGrp: #try magnetic |
|---|
| 202 | MSpGrp = blk.get("_space_group.magn_name_BNS",'') |
|---|
| 203 | if not MSpGrp: |
|---|
| 204 | MSpGrp = blk.get("_space_group_magn.name_BNS",'') |
|---|
| 205 | if not MSpGrp: |
|---|
| 206 | msg = 'No recognizable space group name was found in the CIF.' |
|---|
| 207 | self.errors = msg |
|---|
| 208 | self.warnings += '\n'+msg |
|---|
| 209 | return False |
|---|
| 210 | SpGrp = blk.get('_parent_space_group.name_H-M_alt') |
|---|
| 211 | if not SpGrp: |
|---|
| 212 | SpGrp = blk.get('_parent_space_group.name_H-M') |
|---|
| 213 | # SpGrp = MSpGrp.replace("'",'') |
|---|
| 214 | SpGrp = SpGrp[:2]+SpGrp[2:].replace('_','') #get rid of screw '_' |
|---|
| 215 | if '_' in SpGrp[1]: SpGrp = SpGrp.split('_')[0]+SpGrp[3:] |
|---|
| 216 | SpGrp = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 217 | magnetic = True |
|---|
| 218 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 219 | self.MPhase['General']['AtomPtrs'] = [3,1,10,12] |
|---|
| 220 | if not SpGrp: |
|---|
| 221 | print (MSpGrp) |
|---|
| 222 | self.warnings += 'No space group name was found in the CIF.' |
|---|
| 223 | return False |
|---|
| 224 | else: |
|---|
| 225 | SpGrp = SpGrp.replace('_','') |
|---|
| 226 | self.Phase['General']['Type'] = 'nuclear' |
|---|
| 227 | #process space group symbol |
|---|
| 228 | E,SGData = G2spc.SpcGroup(SpGrp) |
|---|
| 229 | if E and SpGrp: |
|---|
| 230 | SpGrpNorm = G2spc.StandardizeSpcName(SpGrp) |
|---|
| 231 | if SpGrpNorm: |
|---|
| 232 | E,SGData = G2spc.SpcGroup(SpGrpNorm) |
|---|
| 233 | # nope, try the space group "out of the Box" |
|---|
| 234 | if E: |
|---|
| 235 | self.warnings += 'ERROR in space group symbol '+SpGrp |
|---|
| 236 | self.warnings += '\nThe space group has been set to "P 1". ' |
|---|
| 237 | self.warnings += "Change this in phase's General tab." |
|---|
| 238 | self.warnings += '\nAre there spaces separating axial fields?\n\nError msg: ' |
|---|
| 239 | self.warnings += G2spc.SGErrors(E) |
|---|
| 240 | SGData = G2obj.P1SGData # P 1 |
|---|
| 241 | self.Phase['General']['SGData'] = SGData |
|---|
| 242 | if Super: |
|---|
| 243 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 244 | if E: |
|---|
| 245 | self.warnings += 'Invalid super symmetry symbol '+SpGrp+SuperSg |
|---|
| 246 | self.warnings += '\n'+E |
|---|
| 247 | SuperSg = SuperSg[:SuperSg.index(')')+1] |
|---|
| 248 | self.warnings += '\nNew super symmetry symbol '+SpGrp+SuperSg |
|---|
| 249 | E,SSGData = G2spc.SSpcGroup(SGData,SuperSg) |
|---|
| 250 | self.Phase['General']['SSGData'] = SSGData |
|---|
| 251 | if magnetic: |
|---|
| 252 | self.MPhase['General']['SGData'] = SGData |
|---|
| 253 | self.MPhase['General']['SGData']['MagSpGrp'] = MSSpGrp.replace(',','').replace('\\','') |
|---|
| 254 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 255 | |
|---|
| 256 | if magnetic: #replace std operaors with those from cif file - probably not the same! |
|---|
| 257 | SGData['SGFixed'] = True |
|---|
| 258 | SGData['SGOps'] = [] |
|---|
| 259 | SGData['SGCen'] = [] |
|---|
| 260 | if Super: |
|---|
| 261 | SSGData['SSGOps'] = [] |
|---|
| 262 | SSGData['SSGCen'] = [] |
|---|
| 263 | try: |
|---|
| 264 | sgoploop = blk.GetLoop('_space_group_symop_magn_ssg_operation.id') |
|---|
| 265 | sgcenloop = blk.GetLoop('_space_group_symop_magn_ssg_centering.id') |
|---|
| 266 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_ssg_operation.algebraic')[1] |
|---|
| 267 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_ssg_centering.algebraic')[1] |
|---|
| 268 | except KeyError: #old mag cif names |
|---|
| 269 | sgoploop = blk.GetLoop('_space_group_symop.magn_ssg_id') |
|---|
| 270 | sgcenloop = blk.GetLoop('_space_group_symop.magn_ssg_centering_id') |
|---|
| 271 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_ssg_operation_algebraic')[1] |
|---|
| 272 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_ssg_centering_algebraic')[1] |
|---|
| 273 | spnflp = [] |
|---|
| 274 | for op in sgoploop: |
|---|
| 275 | M,T,S = G2spc.MagSSText2MTS(op[opid]) |
|---|
| 276 | SSGData['SSGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 277 | SGData['SGOps'].append([np.array(M,dtype=float)[:3,:3],T[:3]]) |
|---|
| 278 | spnflp.append(S) |
|---|
| 279 | censpn = [] |
|---|
| 280 | for cent in sgcenloop: |
|---|
| 281 | M,C,S = G2spc.MagSSText2MTS(cent[centid]) |
|---|
| 282 | SSGData['SSGCen'].append(C) |
|---|
| 283 | SGData['SGCen'].append(C[:3]) |
|---|
| 284 | censpn += list(np.array(spnflp)*S) |
|---|
| 285 | self.MPhase['General']['SSGData'] = SSGData |
|---|
| 286 | else: |
|---|
| 287 | try: |
|---|
| 288 | sgoploop = blk.GetLoop('_space_group_symop_magn.id') |
|---|
| 289 | sgcenloop = blk.GetLoop('_space_group_symop_magn_centering.id') |
|---|
| 290 | opid = sgoploop.GetItemPosition('_space_group_symop_magn_operation.xyz')[1] |
|---|
| 291 | centid = sgcenloop.GetItemPosition('_space_group_symop_magn_centering.xyz')[1] |
|---|
| 292 | except KeyError: #old mag cif names |
|---|
| 293 | sgoploop = blk.GetLoop('_space_group_symop.magn_id') |
|---|
| 294 | sgcenloop = blk.GetLoop('_space_group_symop.magn_centering_id') |
|---|
| 295 | opid = sgoploop.GetItemPosition('_space_group_symop.magn_operation_xyz')[1] |
|---|
| 296 | centid = sgcenloop.GetItemPosition('_space_group_symop.magn_centering_xyz')[1] |
|---|
| 297 | spnflp = [] |
|---|
| 298 | for op in sgoploop: |
|---|
| 299 | try: |
|---|
| 300 | M,T,S = G2spc.MagText2MTS(op[opid]) |
|---|
| 301 | SGData['SGOps'].append([np.array(M,dtype=float),T]) |
|---|
| 302 | spnflp.append(S) |
|---|
| 303 | except KeyError: |
|---|
| 304 | self.warnings += 'Space group operator '+op[opid]+' is not recognized by GSAS-II' |
|---|
| 305 | return False |
|---|
| 306 | censpn = [] |
|---|
| 307 | for cent in sgcenloop: |
|---|
| 308 | M,C,S = G2spc.MagText2MTS(cent[centid]) |
|---|
| 309 | SGData['SGCen'].append(C) |
|---|
| 310 | censpn += list(np.array(spnflp)*S) |
|---|
| 311 | self.MPhase['General']['SGData'] = SGData |
|---|
| 312 | self.MPhase['General']['SGData']['SpnFlp'] = censpn |
|---|
| 313 | G2spc.GenMagOps(SGData) #set magMom |
|---|
| 314 | self.MPhase['General']['SGData']['MagSpGrp'] = MSpGrp |
|---|
| 315 | MagPtGp = blk.get('_space_group.magn_point_group') |
|---|
| 316 | if not MagPtGp: |
|---|
| 317 | MagPtGp = blk.get('_space_group.magn_point_group_name') |
|---|
| 318 | if not MagPtGp: |
|---|
| 319 | MagPtGp = blk.get('_space_group_magn.point_group_name') |
|---|
| 320 | self.MPhase['General']['SGData']['MagPtGp'] = MagPtGp |
|---|
| 321 | |
|---|
| 322 | # cell parameters |
|---|
| 323 | cell = [] |
|---|
| 324 | for lbl in cellitems: |
|---|
| 325 | cell.append(cif.get_number_with_esd(blk[lbl])[0]) |
|---|
| 326 | Volume = G2lat.calc_V(G2lat.cell2A(cell)) |
|---|
| 327 | self.Phase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 328 | if magnetic: |
|---|
| 329 | self.MPhase['General']['Cell'] = [False,]+cell+[Volume,] |
|---|
| 330 | if Super: |
|---|
| 331 | waveloop = blk.GetLoop('_cell_wave_vector_seq_id') |
|---|
| 332 | waveDict = dict(waveloop.items()) |
|---|
| 333 | SuperVec = [[cif.get_number_with_esd(waveDict['_cell_wave_vector_x'][0].replace('?','0'))[0], |
|---|
| 334 | cif.get_number_with_esd(waveDict['_cell_wave_vector_y'][0].replace('?','0'))[0], |
|---|
| 335 | cif.get_number_with_esd(waveDict['_cell_wave_vector_z'][0].replace('?','0'))[0]],False,4] |
|---|
| 336 | |
|---|
| 337 | # read in atoms |
|---|
| 338 | self.errors = 'Error during reading of atoms' |
|---|
| 339 | atomlbllist = [] # table to look up atom IDs |
|---|
| 340 | atomloop = blk.GetLoop('_atom_site_label') |
|---|
| 341 | atomkeys = [i.lower() for i in atomloop.keys()] |
|---|
| 342 | if not blk.get('_atom_site_type_symbol'): |
|---|
| 343 | self.isodistort_warnings += '\natom types are missing. \n Check & revise atom types as needed' |
|---|
| 344 | if magnetic: |
|---|
| 345 | try: |
|---|
| 346 | magmoment = '_atom_site_moment.label' |
|---|
| 347 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 348 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 349 | magatomlabels = blk.get(magmoment) |
|---|
| 350 | G2MagDict = {'_atom_site_moment.label': 0, |
|---|
| 351 | '_atom_site_moment.crystalaxis_x':7, |
|---|
| 352 | '_atom_site_moment.crystalaxis_y':8, |
|---|
| 353 | '_atom_site_moment.crystalaxis_z':9} |
|---|
| 354 | except KeyError: |
|---|
| 355 | magmoment = '_atom_site_moment_label' |
|---|
| 356 | magatomloop = blk.GetLoop(magmoment) |
|---|
| 357 | magatomkeys = [i.lower() for i in magatomloop.keys()] |
|---|
| 358 | magatomlabels = blk.get(magmoment) |
|---|
| 359 | G2MagDict = {'_atom_site_moment_label': 0, |
|---|
| 360 | '_atom_site_moment_crystalaxis_x':7, |
|---|
| 361 | '_atom_site_moment_crystalaxis_y':8, |
|---|
| 362 | '_atom_site_moment_crystalaxis_z':9} |
|---|
| 363 | |
|---|
| 364 | if blk.get('_atom_site_aniso_label'): |
|---|
| 365 | anisoloop = blk.GetLoop('_atom_site_aniso_label') |
|---|
| 366 | anisokeys = [i.lower() for i in anisoloop.keys()] |
|---|
| 367 | anisolabels = blk.get('_atom_site_aniso_label') |
|---|
| 368 | else: |
|---|
| 369 | anisoloop = None |
|---|
| 370 | anisokeys = [] |
|---|
| 371 | anisolabels = [] |
|---|
| 372 | if Super: |
|---|
| 373 | occFloop = None |
|---|
| 374 | occCloop = None |
|---|
| 375 | occFdict = {} |
|---|
| 376 | occCdict = {} |
|---|
| 377 | displSloop = None |
|---|
| 378 | displFloop = None |
|---|
| 379 | MagFloop = None |
|---|
| 380 | displSdict = {} |
|---|
| 381 | displFdict = {} |
|---|
| 382 | MagFdict = {} |
|---|
| 383 | UijFloop = None |
|---|
| 384 | UijFdict = {} |
|---|
| 385 | #occupancy modulation |
|---|
| 386 | if blk.get('_atom_site_occ_Fourier_atom_site_label'): |
|---|
| 387 | occFloop = blk.GetLoop('_atom_site_occ_Fourier_atom_site_label') |
|---|
| 388 | occFdict = dict(occFloop.items()) |
|---|
| 389 | if blk.get('_atom_site_occ_special_func_atom_site_label'): #Crenel (i.e. Block Wave) occ |
|---|
| 390 | occCloop = blk.GetLoop('_atom_site_occ_special_func_atom_site_label') |
|---|
| 391 | occCdict = dict(occCloop.items()) |
|---|
| 392 | #position modulation |
|---|
| 393 | if blk.get('_atom_site_displace_Fourier_atom_site_label'): |
|---|
| 394 | displFloop = blk.GetLoop('_atom_site_displace_Fourier_atom_site_label') |
|---|
| 395 | displFdict = dict(displFloop.items()) |
|---|
| 396 | if blk.get('_atom_site_displace_special_func_atom_site_label'): #sawtooth |
|---|
| 397 | displSloop = blk.GetLoop('_atom_site_displace_special_func_atom_site_label') |
|---|
| 398 | displSdict = dict(displSloop.items()) |
|---|
| 399 | #U modulation |
|---|
| 400 | if blk.get('_atom_site_U_Fourier_atom_site_label'): |
|---|
| 401 | UijFloop = blk.GetLoop('_atom_site_U_Fourier_atom_site_label') |
|---|
| 402 | UijFdict = dict(UijFloop.items()) |
|---|
| 403 | #Mag moment modulation |
|---|
| 404 | if blk.get('_atom_site_moment_Fourier_atom_site_label'): |
|---|
| 405 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier_atom_site_label') |
|---|
| 406 | MagFdict = dict(MagFloop.items()) |
|---|
| 407 | Mnames = ['_atom_site_moment_fourier_atom_site_label', |
|---|
| 408 | '_atom_site_moment_fourier_axis','_atom_site_moment_fourier_wave_vector_seq_id', |
|---|
| 409 | '_atom_site_moment_fourier_param_sin','_atom_site_moment_fourier_param_cos'] |
|---|
| 410 | elif blk.get('_atom_site_moment_Fourier.atom_site_label'): |
|---|
| 411 | MagFloop = blk.GetLoop('_atom_site_moment_Fourier.atom_site_label') |
|---|
| 412 | MagFdict = dict(MagFloop.items()) |
|---|
| 413 | Mnames = ['_atom_site_moment_fourier.atom_site_label', |
|---|
| 414 | '_atom_site_moment_fourier.axis','_atom_site_moment_fourier.wave_vector_seq_id', |
|---|
| 415 | '_atom_site_moment_fourier_param.sin','_atom_site_moment_fourier_param.cos'] |
|---|
| 416 | self.Phase['Atoms'] = [] |
|---|
| 417 | if magnetic: |
|---|
| 418 | self.MPhase['Atoms'] = [] |
|---|
| 419 | G2AtomDict = { '_atom_site_type_symbol' : 1, |
|---|
| 420 | '_atom_site_label' : 0, |
|---|
| 421 | '_atom_site_fract_x' : 3, |
|---|
| 422 | '_atom_site_fract_y' : 4, |
|---|
| 423 | '_atom_site_fract_z' : 5, |
|---|
| 424 | '_atom_site_occupancy' : 6, |
|---|
| 425 | '_atom_site_aniso_u_11' : 11, |
|---|
| 426 | '_atom_site_aniso_u_22' : 12, |
|---|
| 427 | '_atom_site_aniso_u_33' : 13, |
|---|
| 428 | '_atom_site_aniso_u_12' : 14, |
|---|
| 429 | '_atom_site_aniso_u_13' : 15, |
|---|
| 430 | '_atom_site_aniso_u_23' : 16, } |
|---|
| 431 | |
|---|
| 432 | ranIdlookup = {} |
|---|
| 433 | for aitem in atomloop: |
|---|
| 434 | atomlist = ['','','',0.,0.,0.,1.0,'',0.,'I',0.01, 0.,0.,0.,0.,0.,0.,] |
|---|
| 435 | for val,key in zip(aitem,atomkeys): |
|---|
| 436 | col = G2AtomDict.get(key,-1) |
|---|
| 437 | if col >= 3: |
|---|
| 438 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 439 | if col >= 11: atomlist[9] = 'A' # if any Aniso term is defined, set flag |
|---|
| 440 | elif col is not None and col != -1: |
|---|
| 441 | atomlist[col] = val |
|---|
| 442 | elif key in ('_atom_site_thermal_displace_type', |
|---|
| 443 | '_atom_site_adp_type'): #Iso or Aniso? |
|---|
| 444 | if val.lower() == 'uani': |
|---|
| 445 | atomlist[9] = 'A' |
|---|
| 446 | elif key == '_atom_site_u_iso_or_equiv': |
|---|
| 447 | uisoval = cif.get_number_with_esd(val)[0] |
|---|
| 448 | if uisoval is not None: |
|---|
| 449 | atomlist[10] = uisoval |
|---|
| 450 | if not atomlist[1] and atomlist[0]: |
|---|
| 451 | typ = atomlist[0].rstrip('0123456789-+') |
|---|
| 452 | if G2elem.CheckElement(typ): |
|---|
| 453 | atomlist[1] = typ |
|---|
| 454 | if not atomlist[1]: |
|---|
| 455 | atomlist[1] = 'Xe' |
|---|
| 456 | self.warnings += ' Atom type '+typ+' not recognized; Xe assumed\n' |
|---|
| 457 | if atomlist[0] in anisolabels: # does this atom have aniso values in separate loop? |
|---|
| 458 | atomlist[9] = 'A' |
|---|
| 459 | for val,key in zip(anisoloop.GetKeyedPacket('_atom_site_aniso_label',atomlist[0]),anisokeys): |
|---|
| 460 | col = G2AtomDict.get(key) |
|---|
| 461 | if col: |
|---|
| 462 | atomlist[col] = cif.get_number_with_esd(val)[0] |
|---|
| 463 | atomlist[7],atomlist[8] = G2spc.SytSym(atomlist[3:6],SGData)[:2] |
|---|
| 464 | atomlist[1] = G2elem.FixValence(atomlist[1]) |
|---|
| 465 | atomlist.append(ran.randint(0,sys.maxsize)) # add a random Id |
|---|
| 466 | self.Phase['Atoms'].append(atomlist) |
|---|
| 467 | ranIdlookup[atomlist[0]] = atomlist[-1] |
|---|
| 468 | if atomlist[0] in atomlbllist: |
|---|
| 469 | self.warnings += ' ERROR: repeated atom label: '+atomlist[0] |
|---|
| 470 | else: |
|---|
| 471 | atomlbllist.append(atomlist[0]) |
|---|
| 472 | |
|---|
| 473 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 474 | matomlist = atomlist[:7]+[0.,0.,0.,]+atomlist[7:] |
|---|
| 475 | for mval,mkey in zip(magatomloop.GetKeyedPacket(magmoment,atomlist[0]),magatomkeys): |
|---|
| 476 | mcol = G2MagDict.get(mkey,-1) |
|---|
| 477 | if mcol > 0: |
|---|
| 478 | matomlist[mcol] = cif.get_number_with_esd(mval)[0] |
|---|
| 479 | self.MPhase['Atoms'].append(matomlist) |
|---|
| 480 | if Super: |
|---|
| 481 | Sfrac = np.zeros((4,2)) |
|---|
| 482 | Sadp = np.zeros((4,12)) |
|---|
| 483 | Spos = np.zeros((4,6)) |
|---|
| 484 | Smag = np.zeros((4,6)) |
|---|
| 485 | nim = -1 |
|---|
| 486 | waveType = 'Fourier' |
|---|
| 487 | if occCdict: |
|---|
| 488 | for i,item in enumerate(occCdict['_atom_site_occ_special_func_atom_site_label']): |
|---|
| 489 | if item == atomlist[0]: |
|---|
| 490 | waveType = 'Crenel' |
|---|
| 491 | val = occCdict['_atom_site_occ_special_func_crenel_c'][i] |
|---|
| 492 | Sfrac[0][0] = cif.get_number_with_esd(val)[0] |
|---|
| 493 | val = occCdict['_atom_site_occ_special_func_crenel_w'][i] |
|---|
| 494 | Sfrac[0][1] = cif.get_number_with_esd(val)[0] |
|---|
| 495 | nim = 1 |
|---|
| 496 | |
|---|
| 497 | if nim >= 0: |
|---|
| 498 | Sfrac = [waveType,]+[[sfrac,False] for sfrac in Sfrac[:nim]] |
|---|
| 499 | else: |
|---|
| 500 | Sfrac = [] |
|---|
| 501 | nim = -1 |
|---|
| 502 | if displFdict: |
|---|
| 503 | for i,item in enumerate(displFdict['_atom_site_displace_fourier_atom_site_label']): |
|---|
| 504 | if item == atomlist[0]: |
|---|
| 505 | waveType = 'Fourier' |
|---|
| 506 | ix = ['x','y','z'].index(displFdict['_atom_site_displace_fourier_axis'][i]) |
|---|
| 507 | im = int(displFdict['_atom_site_displace_fourier_wave_vector_seq_id'][i]) |
|---|
| 508 | if im != nim: |
|---|
| 509 | nim = im |
|---|
| 510 | val = displFdict['_atom_site_displace_fourier_param_sin'][i] |
|---|
| 511 | Spos[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 512 | val = displFdict['_atom_site_displace_fourier_param_cos'][i] |
|---|
| 513 | Spos[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 514 | if nim >= 0: |
|---|
| 515 | Spos = [waveType,]+[[spos,False] for spos in Spos[:nim]] |
|---|
| 516 | else: |
|---|
| 517 | Spos = [] |
|---|
| 518 | nim = -1 |
|---|
| 519 | if UijFdict: |
|---|
| 520 | nim = -1 |
|---|
| 521 | for i,item in enumerate(UijFdict['_atom_site_u_fourier_atom_site_label']): |
|---|
| 522 | if item == atomlist[0]: |
|---|
| 523 | ix = ['U11','U22','U33','U12','U13','U23'].index(UijFdict['_atom_site_u_fourier_tens_elem'][i]) |
|---|
| 524 | im = int(UijFdict['_atom_site_u_fourier_wave_vector_seq_id'][i]) |
|---|
| 525 | if im != nim: |
|---|
| 526 | nim = im |
|---|
| 527 | val = UijFdict['_atom_site_u_fourier_param_sin'][i] |
|---|
| 528 | Sadp[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 529 | val = UijFdict['_atom_site_u_fourier_param_cos'][i] |
|---|
| 530 | Sadp[im-1][ix+6] = cif.get_number_with_esd(val)[0] |
|---|
| 531 | if nim >= 0: |
|---|
| 532 | Sadp = ['Fourier',]+[[sadp,False] for sadp in Sadp[:nim]] |
|---|
| 533 | else: |
|---|
| 534 | Sadp = [] |
|---|
| 535 | nim = -1 |
|---|
| 536 | if MagFdict: |
|---|
| 537 | nim = -1 |
|---|
| 538 | for i,item in enumerate(MagFdict[Mnames[0]]): |
|---|
| 539 | if item == atomlist[0]: |
|---|
| 540 | ix = ['x','y','z'].index(MagFdict[Mnames[1]][i]) |
|---|
| 541 | im = int(MagFdict[Mnames[2]][i]) |
|---|
| 542 | if im != nim: |
|---|
| 543 | nim = im |
|---|
| 544 | val = MagFdict[Mnames[3]][i] |
|---|
| 545 | Smag[im-1][ix] = cif.get_number_with_esd(val)[0] |
|---|
| 546 | val = MagFdict[Mnames[4]][i] |
|---|
| 547 | Smag[im-1][ix+3] = cif.get_number_with_esd(val)[0] |
|---|
| 548 | if nim >= 0: |
|---|
| 549 | Smag = ['Fourier',]+[[smag,False] for smag in Smag[:nim]] |
|---|
| 550 | else: |
|---|
| 551 | Smag = [] |
|---|
| 552 | SSdict = {'SS1':{'Sfrac':Sfrac,'Spos':Spos,'Sadp':Sadp,'Smag':Smag}} |
|---|
| 553 | if magnetic and atomlist[0] in magatomlabels: |
|---|
| 554 | matomlist.append(SSdict) |
|---|
| 555 | atomlist.append(SSdict) |
|---|
| 556 | if len(atomlbllist) != len(self.Phase['Atoms']): |
|---|
| 557 | self.isodistort_warnings += '\nRepeated atom labels prevents ISODISTORT decode' |
|---|
| 558 | for lbl in phasenamefields: # get a name for the phase |
|---|
| 559 | name = blk.get(lbl) |
|---|
| 560 | if name is None: |
|---|
| 561 | continue |
|---|
| 562 | name = name.strip() |
|---|
| 563 | if name == '?' or name == '.': |
|---|
| 564 | continue |
|---|
| 565 | else: |
|---|
| 566 | break |
|---|
| 567 | else: # no name found, use block name for lack of a better choice |
|---|
| 568 | name = blknm |
|---|
| 569 | self.Phase['General']['Name'] = name.strip()[:20] |
|---|
| 570 | self.Phase['General']['Super'] = Super |
|---|
| 571 | if magnetic: |
|---|
| 572 | self.MPhase['General']['Type'] = 'magnetic' |
|---|
| 573 | self.MPhase['General']['Name'] = name.strip()[:20]+' mag' |
|---|
| 574 | self.MPhase['General']['Super'] = Super |
|---|
| 575 | if Super: |
|---|
| 576 | if self.MPhase['General']['SGData']['SGGray']: |
|---|
| 577 | SuperSg += 's' |
|---|
| 578 | self.MPhase['General']['Modulated'] = True |
|---|
| 579 | self.MPhase['General']['SuperVec'] = SuperVec |
|---|
| 580 | self.MPhase['General']['SuperSg'] = SuperSg |
|---|
| 581 | else: |
|---|
| 582 | self.MPhase = None |
|---|
| 583 | if Super: |
|---|
| 584 | self.Phase['General']['Modulated'] = True |
|---|
| 585 | self.Phase['General']['SuperVec'] = SuperVec |
|---|
| 586 | self.Phase['General']['SuperSg'] = SuperSg |
|---|
| 587 | if not self.Phase['General']['SGData']['SGFixed']: |
|---|
| 588 | self.Phase['General']['SSGData'] = G2spc.SSpcGroup(SGData,SuperSg)[1] |
|---|
| 589 | if not self.isodistort_warnings: |
|---|
| 590 | if blk.get('_iso_displacivemode_label') or blk.get('_iso_occupancymode_label'): |
|---|
| 591 | self.errors = "Error while processing ISODISTORT constraints" |
|---|
| 592 | self.ISODISTORT_proc(blk,atomlbllist,ranIdlookup) |
|---|
| 593 | else: |
|---|
| 594 | self.warnings += self.isodistort_warnings |
|---|
| 595 | returnstat = True |
|---|
| 596 | return returnstat |
|---|
| 597 | |
|---|
| 598 | def ISODISTORT_proc(self,blk,atomlbllist,ranIdlookup): |
|---|
| 599 | 'Process ISODISTORT items to create constraints etc.' |
|---|
| 600 | varLookup = {'dx':'dAx','dy':'dAy','dz':'dAz','do':'Afrac'} |
|---|
| 601 | 'Maps ISODISTORT parm names to GSAS-II names' |
|---|
| 602 | #---------------------------------------------------------------------- |
|---|
| 603 | # read in the ISODISTORT displacement modes |
|---|
| 604 | #---------------------------------------------------------------------- |
|---|
| 605 | self.Constraints = [] |
|---|
| 606 | explaination = {} |
|---|
| 607 | if blk.get('_iso_displacivemode_label'): |
|---|
| 608 | modelist = [] |
|---|
| 609 | shortmodelist = [] |
|---|
| 610 | for lbl in blk.get('_iso_displacivemode_label'): |
|---|
| 611 | modelist.append(lbl) |
|---|
| 612 | # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB |
|---|
| 613 | # where SSSSS is the parent spacegroup, [x,y,z] is a location |
|---|
| 614 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl) |
|---|
| 615 | # this extracts the AAAAA and BBBBB parts of the string |
|---|
| 616 | if regexp: |
|---|
| 617 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 618 | G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list |
|---|
| 619 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 620 | coordVarLbl = [] |
|---|
| 621 | G2varLbl = [] |
|---|
| 622 | G2varObj = [] |
|---|
| 623 | error = False |
|---|
| 624 | for lbl in blk.get('_iso_deltacoordinate_label'): |
|---|
| 625 | coordVarLbl.append(lbl) |
|---|
| 626 | if '_' in lbl: |
|---|
| 627 | albl = lbl[:lbl.rfind('_')] |
|---|
| 628 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 629 | else: |
|---|
| 630 | self.warnings += ' ERROR: _iso_deltacoordinate_label not parsed: '+lbl |
|---|
| 631 | error = True |
|---|
| 632 | continue |
|---|
| 633 | if albl not in atomlbllist: |
|---|
| 634 | self.warnings += ' ERROR: _iso_deltacoordinate_label atom not found: '+lbl |
|---|
| 635 | error = True |
|---|
| 636 | continue |
|---|
| 637 | else: |
|---|
| 638 | anum = atomlbllist.index(albl) |
|---|
| 639 | var = varLookup.get(vlbl) |
|---|
| 640 | if not var: |
|---|
| 641 | self.warnings += ' ERROR: _iso_deltacoordinate_label variable not found: '+lbl |
|---|
| 642 | error = True |
|---|
| 643 | continue |
|---|
| 644 | G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID |
|---|
| 645 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 646 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 647 | )) |
|---|
| 648 | if error: |
|---|
| 649 | raise Exception("Error decoding variable labels") |
|---|
| 650 | |
|---|
| 651 | if len(G2varObj) != len(modelist): |
|---|
| 652 | print ("non-square input") |
|---|
| 653 | raise Exception("Rank of _iso_displacivemode != _iso_deltacoordinate") |
|---|
| 654 | |
|---|
| 655 | error = False |
|---|
| 656 | ParentCoordinates = {} |
|---|
| 657 | for lbl,exp in zip( |
|---|
| 658 | blk.get('_iso_coordinate_label'), |
|---|
| 659 | blk.get('_iso_coordinate_formula'), |
|---|
| 660 | ): |
|---|
| 661 | if '_' in lbl: |
|---|
| 662 | albl = lbl[:lbl.rfind('_')] |
|---|
| 663 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 664 | else: |
|---|
| 665 | self.warnings += ' ERROR: _iso_coordinate_label not parsed: '+lbl |
|---|
| 666 | error = True |
|---|
| 667 | continue |
|---|
| 668 | if vlbl not in 'xyz' or len(vlbl) != 1: |
|---|
| 669 | self.warnings += ' ERROR: _iso_coordinate_label coordinate not parsed: '+lbl |
|---|
| 670 | error = True |
|---|
| 671 | continue |
|---|
| 672 | i = 'xyz'.index(vlbl) |
|---|
| 673 | if not ParentCoordinates.get(albl): |
|---|
| 674 | ParentCoordinates[albl] = [None,None,None] |
|---|
| 675 | if '+' in exp: |
|---|
| 676 | val = exp.split('+')[0].strip() |
|---|
| 677 | val = G2p3.FormulaEval(val) |
|---|
| 678 | if val is None: |
|---|
| 679 | self.warnings += ' ERROR: _iso_coordinate_formula coordinate not interpreted: '+lbl |
|---|
| 680 | error = True |
|---|
| 681 | continue |
|---|
| 682 | ParentCoordinates[albl][i] = val |
|---|
| 683 | else: |
|---|
| 684 | ParentCoordinates[albl][i] = G2p3.FormulaEval(exp) |
|---|
| 685 | if error: |
|---|
| 686 | print (self.warnings) |
|---|
| 687 | raise Exception("Error decoding variable labels") |
|---|
| 688 | # get mapping of modes to atomic coordinate displacements |
|---|
| 689 | displacivemodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 690 | for row,col,val in zip( |
|---|
| 691 | blk.get('_iso_displacivemodematrix_row'), |
|---|
| 692 | blk.get('_iso_displacivemodematrix_col'), |
|---|
| 693 | blk.get('_iso_displacivemodematrix_value'),): |
|---|
| 694 | displacivemodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 695 | # Invert to get mapping of atom displacements to modes |
|---|
| 696 | displacivemodeInvmatrix = np.linalg.inv(displacivemodematrix) |
|---|
| 697 | # create the constraints |
|---|
| 698 | for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 699 | constraint = [] |
|---|
| 700 | for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 701 | if k == 0: continue |
|---|
| 702 | constraint.append([k,G2varObj[j]]) |
|---|
| 703 | constraint += [shortmodelist[i],False,'f'] |
|---|
| 704 | self.Constraints.append(constraint) |
|---|
| 705 | #---------------------------------------------------------------------- |
|---|
| 706 | # save the ISODISTORT info for "mode analysis" |
|---|
| 707 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 708 | self.Phase['ISODISTORT'].update({ |
|---|
| 709 | 'IsoModeList' : modelist, |
|---|
| 710 | 'G2ModeList' : shortmodelist, |
|---|
| 711 | 'IsoVarList' : coordVarLbl, |
|---|
| 712 | 'G2VarList' : G2varObj, |
|---|
| 713 | 'ParentStructure' : ParentCoordinates, |
|---|
| 714 | 'Var2ModeMatrix' : displacivemodeInvmatrix, |
|---|
| 715 | 'Mode2VarMatrix' : displacivemodematrix, |
|---|
| 716 | }) |
|---|
| 717 | # make explaination dictionary |
|---|
| 718 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 719 | explaination[shortmode] = "ISODISTORT full name "+str(mode) |
|---|
| 720 | #---------------------------------------------------------------------- |
|---|
| 721 | # now read in the ISODISTORT occupancy modes |
|---|
| 722 | #---------------------------------------------------------------------- |
|---|
| 723 | if blk.get('_iso_occupancymode_label'): |
|---|
| 724 | modelist = [] |
|---|
| 725 | shortmodelist = [] |
|---|
| 726 | for lbl in blk.get('_iso_occupancymode_label'): |
|---|
| 727 | modelist.append(lbl) |
|---|
| 728 | # assume lbl is of form SSSSS[x,y,z]AAAA(a,b,...)BBBBB |
|---|
| 729 | # where SSSSS is the parent spacegroup, [x,y,z] is a location |
|---|
| 730 | regexp = re.match(r'.*?\[.*?\](.*?)\(.*?\)(.*)',lbl) |
|---|
| 731 | # this extracts the AAAAA and BBBBB parts of the string |
|---|
| 732 | if regexp: |
|---|
| 733 | lbl = regexp.expand(r'\1\2') # parse succeeded, make a short version |
|---|
| 734 | lbl = lbl.replace('order','o') |
|---|
| 735 | G2obj.MakeUniqueLabel(lbl,shortmodelist) # make unique and add to list |
|---|
| 736 | # read in the coordinate offset variables names and map them to G2 names/objects |
|---|
| 737 | occVarLbl = [] |
|---|
| 738 | G2varLbl = [] |
|---|
| 739 | G2varObj = [] |
|---|
| 740 | error = False |
|---|
| 741 | for lbl in blk.get('_iso_deltaoccupancy_label'): |
|---|
| 742 | occVarLbl.append(lbl) |
|---|
| 743 | if '_' in lbl: |
|---|
| 744 | albl = lbl[:lbl.rfind('_')] |
|---|
| 745 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 746 | else: |
|---|
| 747 | self.warnings += ' ERROR: _iso_deltaoccupancy_label not parsed: '+lbl |
|---|
| 748 | error = True |
|---|
| 749 | continue |
|---|
| 750 | if albl not in atomlbllist: |
|---|
| 751 | self.warnings += ' ERROR: _iso_deltaoccupancy_label atom not found: '+lbl |
|---|
| 752 | error = True |
|---|
| 753 | continue |
|---|
| 754 | else: |
|---|
| 755 | anum = atomlbllist.index(albl) |
|---|
| 756 | var = varLookup.get(vlbl) |
|---|
| 757 | if not var: |
|---|
| 758 | self.warnings += ' ERROR: _iso_deltaoccupancy_label variable not found: '+lbl |
|---|
| 759 | error = True |
|---|
| 760 | continue |
|---|
| 761 | G2varLbl.append('::'+var+':'+str(anum)) # variable name, less phase ID |
|---|
| 762 | G2varObj.append(G2obj.G2VarObj( |
|---|
| 763 | (self.Phase['ranId'],None,var,ranIdlookup[albl]) |
|---|
| 764 | )) |
|---|
| 765 | if error: |
|---|
| 766 | raise Exception("Error decoding variable labels") |
|---|
| 767 | |
|---|
| 768 | if len(G2varObj) != len(modelist): |
|---|
| 769 | print ("non-square input") |
|---|
| 770 | raise Exception("Rank of _iso_occupancymode != _iso_deltaoccupancy") |
|---|
| 771 | |
|---|
| 772 | error = False |
|---|
| 773 | ParentCoordinates = {} |
|---|
| 774 | for lbl,exp in zip( |
|---|
| 775 | blk.get('_iso_occupancy_label'), |
|---|
| 776 | blk.get('_iso_occupancy_formula'), |
|---|
| 777 | ): |
|---|
| 778 | if '_' in lbl: |
|---|
| 779 | albl = lbl[:lbl.rfind('_')] |
|---|
| 780 | vlbl = lbl[lbl.rfind('_')+1:] |
|---|
| 781 | else: |
|---|
| 782 | self.warnings += ' ERROR: _iso_occupancy_label not parsed: '+lbl |
|---|
| 783 | error = True |
|---|
| 784 | continue |
|---|
| 785 | if vlbl != 'occ': |
|---|
| 786 | self.warnings += ' ERROR: _iso_occupancy_label coordinate not parsed: '+lbl |
|---|
| 787 | error = True |
|---|
| 788 | continue |
|---|
| 789 | if '+' in exp: |
|---|
| 790 | val = exp.split('+')[0].strip() |
|---|
| 791 | val = G2p3.FormulaEval(val) |
|---|
| 792 | if val is None: |
|---|
| 793 | self.warnings += ' ERROR: _iso_occupancy_formula coordinate not interpreted: '+lbl |
|---|
| 794 | error = True |
|---|
| 795 | continue |
|---|
| 796 | ParentCoordinates[albl] = val |
|---|
| 797 | if error: |
|---|
| 798 | raise Exception("Error decoding occupancy labels") |
|---|
| 799 | # get mapping of modes to atomic coordinate displacements |
|---|
| 800 | occupancymodematrix = np.zeros((len(G2varObj),len(G2varObj))) |
|---|
| 801 | for row,col,val in zip( |
|---|
| 802 | blk.get('_iso_occupancymodematrix_row'), |
|---|
| 803 | blk.get('_iso_occupancymodematrix_col'), |
|---|
| 804 | blk.get('_iso_occupancymodematrix_value'),): |
|---|
| 805 | occupancymodematrix[int(row)-1,int(col)-1] = float(val) |
|---|
| 806 | # Invert to get mapping of atom displacements to modes |
|---|
| 807 | occupancymodeInvmatrix = np.linalg.inv(occupancymodematrix) |
|---|
| 808 | # create the constraints |
|---|
| 809 | for i,row in enumerate(occupancymodeInvmatrix): |
|---|
| 810 | constraint = [] |
|---|
| 811 | for j,(lbl,k) in enumerate(zip(occVarLbl,row)): |
|---|
| 812 | if k == 0: continue |
|---|
| 813 | constraint.append([k,G2varObj[j]]) |
|---|
| 814 | constraint += [shortmodelist[i],False,'f'] |
|---|
| 815 | self.Constraints.append(constraint) |
|---|
| 816 | #---------------------------------------------------------------------- |
|---|
| 817 | # save the ISODISTORT info for "mode analysis" |
|---|
| 818 | if 'ISODISTORT' not in self.Phase: self.Phase['ISODISTORT'] = {} |
|---|
| 819 | self.Phase['ISODISTORT'].update({ |
|---|
| 820 | 'OccModeList' : modelist, |
|---|
| 821 | 'G2OccModeList' : shortmodelist, |
|---|
| 822 | 'OccVarList' : occVarLbl, |
|---|
| 823 | 'G2OccVarList' : G2varObj, |
|---|
| 824 | 'BaseOcc' : ParentCoordinates, |
|---|
| 825 | 'Var2OccMatrix' : occupancymodeInvmatrix, |
|---|
| 826 | 'Occ2VarMatrix' : occupancymodematrix, |
|---|
| 827 | }) |
|---|
| 828 | # make explaination dictionary |
|---|
| 829 | for mode,shortmode in zip(modelist,shortmodelist): |
|---|
| 830 | explaination[shortmode] = "ISODISTORT full name "+str(mode) |
|---|
| 831 | #---------------------------------------------------------------------- |
|---|
| 832 | # done with read |
|---|
| 833 | #---------------------------------------------------------------------- |
|---|
| 834 | if explaination: self.Constraints.append(explaination) |
|---|
| 835 | |
|---|
| 836 | # # debug: show the mode var to mode relations |
|---|
| 837 | # for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 838 | # l = '' |
|---|
| 839 | # for j,(lbl,k) in enumerate(zip(coordVarLbl,row)): |
|---|
| 840 | # if k == 0: continue |
|---|
| 841 | # if l: l += ' + ' |
|---|
| 842 | # #l += lbl+' * '+str(k) |
|---|
| 843 | # l += G2varLbl[j]+' * '+str(k) |
|---|
| 844 | # print str(i) + ': '+shortmodelist[i]+' = '+l |
|---|
| 845 | # print 70*'=' |
|---|
| 846 | |
|---|
| 847 | # # debug: Get the ISODISTORT offset values |
|---|
| 848 | # coordVarDelta = {} |
|---|
| 849 | # for lbl,val in zip( |
|---|
| 850 | # blk.get('_iso_deltacoordinate_label'), |
|---|
| 851 | # blk.get('_iso_deltacoordinate_value'),): |
|---|
| 852 | # coordVarDelta[lbl] = float(val) |
|---|
| 853 | # modeVarDelta = {} |
|---|
| 854 | # for lbl,val in zip( |
|---|
| 855 | # blk.get('_iso_displacivemode_label'), |
|---|
| 856 | # blk.get('_iso_displacivemode_value'),): |
|---|
| 857 | # modeVarDelta[lbl] = cif.get_number_with_esd(val)[0] |
|---|
| 858 | |
|---|
| 859 | # print 70*'=' |
|---|
| 860 | # # compute the mode values from the reported coordinate deltas |
|---|
| 861 | # for i,row in enumerate(displacivemodeInvmatrix): |
|---|
| 862 | # l = '' |
|---|
| 863 | # sl = '' |
|---|
| 864 | # s = 0. |
|---|
| 865 | # for lbl,k in zip(coordVarLbl,row): |
|---|
| 866 | # if k == 0: continue |
|---|
| 867 | # if l: l += ' + ' |
|---|
| 868 | # l += lbl+' * '+str(k) |
|---|
| 869 | # if sl: sl += ' + ' |
|---|
| 870 | # sl += str(coordVarDelta[lbl])+' * '+str(k) |
|---|
| 871 | # s += coordVarDelta[lbl] * k |
|---|
| 872 | # print 'a'+str(i)+' = '+l |
|---|
| 873 | # print '\t= '+sl |
|---|
| 874 | # print modelist[i],shortmodelist[i],modeVarDelta[modelist[i]],s |
|---|
| 875 | # print |
|---|
| 876 | |
|---|
| 877 | # print 70*'=' |
|---|
| 878 | # # compute the coordinate displacements from the reported mode values |
|---|
| 879 | # for i,lbl,row in zip(range(len(coordVarLbl)),coordVarLbl,displacivemodematrix): |
|---|
| 880 | # l = '' |
|---|
| 881 | # sl = '' |
|---|
| 882 | # s = 0.0 |
|---|
| 883 | # for j,k in enumerate(row): |
|---|
| 884 | # if k == 0: continue |
|---|
| 885 | # if l: l += ' + ' |
|---|
| 886 | # l += 'a'+str(j+1)+' * '+str(k) |
|---|
| 887 | # if sl: sl += ' + ' |
|---|
| 888 | # sl += str(shortmodelist[j]) +' = '+ str(modeVarDelta[modelist[j]]) + ' * '+str(k) |
|---|
| 889 | # s += modeVarDelta[modelist[j]] * k |
|---|
| 890 | # print lbl+' = '+l |
|---|
| 891 | # print '\t= '+sl |
|---|
| 892 | # print lbl,G2varLbl[i],coordVarDelta[lbl],s |
|---|
| 893 | # print |
|---|
| 894 | |
|---|
| 895 | # determine the coordinate delta values from deviations from the parent structure |
|---|
| 896 | # for atmline in self.Phase['Atoms']: |
|---|
| 897 | # lbl = atmline[0] |
|---|
| 898 | # x,y,z = atmline[3:6] |
|---|
| 899 | # if lbl not in ParentCoordinates: |
|---|
| 900 | # print lbl,x,y,z |
|---|
| 901 | # continue |
|---|
| 902 | # px,py,pz = ParentCoordinates[lbl] |
|---|
| 903 | # print lbl,x,y,z,x-px,y-py,z-pz |
|---|
